% 4th order Runge-Kutta method
function [x,y]=RungeKuttaFn(x0,y0,h,N,f)

x=x0+(0:N)'*h;
y=zeros(N+1,1);
y(1)=y0;
for k=1:N
    k1=f(x(k),y(k));
    k2=f(x(k)+h/2,y(k)+k1*h/2);
    k3=f(x(k)+h/2,y(k)+k2*h/2);
    k4=f(x(k+1),y(k)+k3*h);
    y(k+1)=y(k)+h/6*(k1+2*k2+2*k3+k4);
end
