function [y time]=rk2(f,y0,t0,dt,nt)
time(1)=t0;
y(:,1)=y0;
nstage=2;
a=[0 0;1 0];
c=[0; 1];
b=[0.5 0.5];
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