function [t,x]=rk23(f,x0,t0,tf,rtol,atol)
   t=zeros(1,10000);
   t(1)=t0; 
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   h=1e-6;
   k=1;
   while t(k)<tf
       if t(k)+h>tf then 
           h=tf-t(k);
       end;
       
       k1=f(x(:,k),t(k));
       k2=f(x(:,k)+h*k1,t(k)+h);
       k3=f(x(:,k)+h*k1/4+h*k2/4,t(k)+h/2);
       xheun=x(:,k)+h*(k1/2+k2/2);
       xrk3=x(:,k)+h*(k1/6+k2/6+2*k3/3);
       err=norm(xheun-xrk3);
       err0=rtol*norm(xrk3)+atol;
       if err<err0 then
           x(:,k+1)=xrk3;
           t(k+1)=t(k)+h;
           k=k+1;
       end
       h=h*(err0/err)^(1/3)*0.8;           
   end
   t=t(1:k);
   x=x(:,1:k);
endfunction
