function [ynew] = rk4(y,t,h,f)
   k1 = f(y,t);
   k2 = f(y + 0.5*h*k1, t + 0.5*h);
   k3 = f(y + 0.5*h*k2, t + 0.5*h);
   k4 = f(y + 1.0*h*k3, t + 1.0*h);
   
   ynew = y + 1.0/6.0*h*(k1+2*k2+2*k3+k4)
end

