function [tout, yout] = solveRK4( yp, tspan, y0, h )
% function [tout,yout]=solveRK4( yp, tspan, y0, deltat ) 
% solves y' = yp using RK4 with step h.
%
% yp, tspan, y0 are as usual in Matlab ODE syntax, tho' tspan is
% assumed to [ t0 tfinal ] since output is spaced at h

%Get initial and final times:
[m n] = size( tspan );
t = tspan(1, 1);
tfinal = tspan(m, n);

%Initialize the differential equation
y = y0(:);
tout = t;
yout = y.';

while (t < tfinal)

   if t + h > tfinal, h = tfinal - t; end

   % Compute the slope
   s1 = feval(yp, t, y); s1 = s1(:);
   s2 = feval(yp, t + h/2, y + h*s1/2); s2 = s2(:);
   s3 = feval(yp, t + h/2, y + h*s2/2); s3 = s3(:);
   s4 = feval(yp, t + h, y + h*s3); s4 = s4(:);

   % Update
   t = t + h;
   n = n + 1;
   y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
   tout = [tout; t];
   yout = [yout; y.'];

end;   % while