function [time,y]=rk4(time,y,dt);
k1=dt*f(time,y);
k2=dt*f(time+dt/2,y+k1/2);
k3=dt*f(time+dt/2,y+k2/2);
k4=dt*f(time+dt,y+k3);
y=y+(k1+2*k2+2*k3+k4)/6;
time=time+dt;
