function [X,Y] = rk(f,span,y,n)
h = (span(2)-span(1))/n;     % step size
x = span(1);                 % initial x
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