function [t,x,r]=rk23ev(fun,x0,d,t0,tf,reltol)
    [f,zc,hpos,hneg]=fun();
    abstol=1e-9;
    hmin=1e-9*(tf-t0);
    hmax=(tf-t0)/50;   
    h=hmin;
    N=2;
    r=0;
    t=zeros(1,100000);
    x=zeros(length(x0),length(t));
    t(1)=t0;
    x(:,1)=x0;
    k=1;
    zcold=zc(x0,d,t0);//funcion de cruce original
    while t(k)<tf
        [x1,x2]=rk23ev_step(f,x(:,k),d,t(k),h);
        e=norm(x1-x2)+1e-12;
        tol=reltol*norm(x1)+abstol;
        if e<tol then
            //verificar cruce por cero
            zcnew=zc(x2,d,t(k)+h);
            if min(zcnew.*zcold)<0 then
                //discontinuidad
                //buscar discontinuidad. Usamos bisección
                tlow=t(k);
                tmax=tlow+h;
                xmax=x2;
                while tmax-tlow>10*hmin
                    tmean=(tlow+tmax)/2;
                    xnew=heunev_step(f,x(:,k),d,t(k),tmean-t(k));
                    zcnew=zc(xnew,d,tmean);
                    if min(zcnew.*zcold)<0 then
                        //tmean está a la derecha de la discontinuidad
                        tmax=tmean;
                        xmax=xnew;
                        lastchange=zcnew.*zcold;
                    else
                        tlow=tmean;
                    end
                end
                x(:,k+1)=xmax;//me quedo con el ultimo punto a la derecha
                t(k+1)=tmax;
                k=k+1;
                //ejecutar los handlers
                for j=1:length(zcold)
                    if zcold(j)<0 & lastchange(j)<0 then
                        [xnew,d]=hpos(x(:,k),d,t(k),j);     
                    elseif zcold(j)>0 & lastchange(j)<0 then
                        [xnew,d]=hneg(x(:,k),d,t(k),j);     
                    end
                end
                zcold=zc(xnew,d,t(k));
                if xnew<>x(:,k) then              
                    x(:,k+1)=xnew;
                    t(k+1)=t(k);
                    k=k+1;
                end
                h=hmin;
            else
                x(:,k+1)=x2;
                t(k+1)=t(k)+h;
                k=k+1;
            end
        else
            r=r+1;    
        end
        h=0.8*h*(tol/e)^(1/(N+1));
        if h<hmin then 
            h=hmin;
        end
        if t(k)+h>tf then 
            h=tf-t(k);
        end
        if h>hmax then 
            h=hmax;
        end
        
    end
    t=t(1:k);
    x=x(:,1:k);
endfunction

function [x1,x2]=rk23ev_step(f,xk,d,tk,h)
    k1=f(xk,d,tk);
    k2=f(xk+h*k1,d,tk+h);
    k3=f(xk+h*k1/4+h*k2/4,d,tk+h/2);
    x1=xk+h*(k1/2+k2/2);
    x2=xk+h*(k1/6+k2/6+k3*2/3);
endfunction

function xk1=heunev_step(f,xk,d,tk,h)
    k1=f(xk,d,tk);
    k2=f(xk+h*k1,d,tk+h);
    xk1=xk+h/2*(k1+k2);
endfunction    
            
    