function [X,Y]=rk(x,y,x1,n)

h=(x1-x)/n;                 % step size
X=x;                        % initial x
Y=y;                        % initial y

for i=1:n                   % begin loop
   
   k1=f(x,y);               % first slope
   k2=f(x+h/2,y+h*k1/2);    % second slope
   k3=f(x+h/2,y+h*k2/2);    % third slope
   k4=f(x+h,y+h*k3);        % fourth slope
   
   k=(k1+2*k2+2*k3+k4)/6;   % average slope
   
   x=x+h;                   % new x
   
   y=y+h*k;                 % new y
   
   X=[X,x];                 % update x-column
   Y=[Y,y];                 % update y-column
end                         % end loop