function [y time kk zz]=rk4st(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];
kk=0;
nd=size(y0);
zz=zeros(9*nd,1);
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;
 if ((y(1,n+1)-1)*(y(1,n)-1)<0)
  kk=kk+1;
  zz((kk-1)*nd+1:kk*nd)=y(:,n+1);
 endif
 if (kk==9) break;
 endif
 if (y(2,n+1)>200) break;
 endif
endfor
end