function [x,y]=Runge_Kutta4(ydot_fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    k1=h*feval(ydot_fun,x(n),y(n));
    k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1);
    k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2);
    k4=h*feval(ydot_fun,x(n)+h,y(n)+k3);
    y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end