function [y time]=rk4(f,y0,t0,dt,nt)
time(1)=t0;
y(:,1)=y0;
nstage=4;
a=[0 0 0 0;1/2 0 0 0;0 1/2 0 0;0 0 1 0];
c=[0; 1/2; 1/2; 1];
b=[1/6 2/6 2/6 1/6];
for n=1:nt
 K(:,1)=feval(f,y(:,n),time(n));
 for m=2:nstage
  tt=time(n)+c(m)*dt;
  yy=y(:,n)+dt*K(:,1:m-1)*a(m,1:m-1)';
  K(:,m)=feval(f,yy,tt);
 endfor
 y(:,n+1)=y(:,n)+dt*K(:,1:nstage)*b';
 time(n+1)=time(n)+dt;
endfor