rk4

function y = rk4(a,b,f,n,y0)
   h = (b-a)/n;
   x = a:h:b;
   y = zeros(length(x),1);
   y(1) = y0;

   for i = 2:length(x)
       s1 = h*feval(f,x(i-1),y(i-1));
       s2 = h*feval(f,x(i-1)+h/2,y(i-1)+s1/2);
       s3 = h*feval(f,x(i-1)+h/2,y(i-1)+s2/2);
       s4 = h*feval(f,x(i-1)+h,y(i-1)+s3);
       y(i) = y(i-1) + 1/6*(s1+2*s2+2*s3+s4);
   end

end