%
% rk4 -- a 4th order Runge Kutta Method
%
% [tv,yv,n]=rk4(a,b,y0,h)
%
function [tv,yv,n]=rk4(a,b,y0,h);
n=floor((b-a)/h);
yv(1)=y0;
tv(1)=a;
for ii=2:n+1;
    f1=fun(tv(ii-1),yv(ii-1));
    tempy=yv(ii-1)+h/2*f1;
    f2=fun(tv(ii-1)+h/2,tempy);
    tempy=yv(ii-1)+h/2*f2;
    f3=fun(tv(ii-1)+h/2,tempy);
    tempy=yv(ii-1)+h*f3;
    f4=fun(tv(ii-1)+h,tempy);
    yv(ii)=yv(ii-1)+h*(f1+2*f2+2*f3+f4)/6;
    tv(ii)=tv(ii-1)+h;
end;