%This is rk4.m

function [tout, yout] = rk4(FunFcn, t0, tfinal, y0, ssize)

if nargin < 5, ssize = (tfinal - t0)/100; end
h=ssize;
t = t0; y=y0(:);
tout = t;
yout = y.';

while (t < tfinal)
      if t + h > tfinal, h = tfinal - t; end

      s1 = feval(FunFcn, t, y); s1 = s1(:);
      s2 = feval(FunFcn, t + h/2, y + h*s1/2); s2=s2(:);
      s3 = feval(FunFcn, t + h/2, y + h*s2/2); s3=s3(:);
      s4 = feval(FunFcn, t + h, y + h*s3); s4=s4(:);
      t = t + h;
      y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6;
      tout = [tout; t];
      yout = [yout; y.'];

end;