function [t,x]=rk23ev(sys,x0,d,t0,tf,rtol,atol)
   t=zeros(1,20000);
   t(1)=t0; 
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   hmax=(tf-t0)/500;
   h=1e-6;
   k=1;
   [f,zc,H]=sys();
   zold=zc(x0,d,t0);
   while t(k)<tf
       if t(k)+h>tf then 
           h=tf-t(k);
       end;
       
       k1=f(x(:,k),d,t(k));
       k2=f(x(:,k)+h*k1,d,t(k)+h);
       k3=f(x(:,k)+h*k1/4+h*k2/4,d,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=max(rtol*norm(xrk3),atol);
       if err<err0 then
           z=zc(xrk3,d,t(k)+h);
           if sign(z)==sign(zold) then
               //do nothing
           else
               
               //discontinuity
               dt=h;
               t1=t(k);
               t2=t(k)+h;
               z1=zold;
               z2=z;
               h1=h;
               while dt>h1/10000;
                   dt=dt/2;
                   h=t1+dt-t(k);
                   k1=f(x(:,k),d,t(k));
                   k2=f(x(:,k)+h*k1,d,t(k)+h);
                   k3=f(x(:,k)+h*k1/4+h*k2/4,d,t(k)+h/2);
                   xnew=x(:,k)+h*(k1/6+k2/6+2*k3/3);
                   z=zc(xnew,d,t(k)+h);
                   if sign(z)==sign(z1) then
                       z1=z;
                       t1=t(k)+h;
                   else
                       z2=z;
                       t2=t(k)+h;
                       xrk3=xnew;  
                   end
               end
               h=t2-t(k);
               for i=1:length(zold)
                   if sign(zold(i))<>sign(z2(i)) then
                       [xrk3,d]=H(xrk3,d,t(k)+h,sign(z2(i)),i);
                   end
               end    
               zold=z2;
           end
           x(:,k+1)=xrk3;
           t(k+1)=t(k)+h;
           k=k+1;
       end
       if err>0 then
           h=h*(err0/err)^(1/3)*0.8;   
           if h>hmax then 
               h=hmax;
           end
                   
       else
           h=hmax;
       end
        
   end
   t=t(1:k);
   x=x(:,1:k);
endfunction
