function [x1,x2]=rkf45sev_step(f,xk,tk,d,h)
  k1=f(xk,tk,d);
  k2=f(xk+1/4*h*k1,tk+1/4*h,d);
  k3=f(xk+3/32*h*k1+9/32*h*k2,tk+3/8*h,d); 
  k4=f(xk+1932/2197*h*k1-7200/2197*h*k2+7296/2197*h*k3,tk+12/13*h,d); 
  k5=f(xk+439/216*h*k1-8*h*k2+3680/513*h*k3-845/4104*h*k4,tk+h,d); 
  k6=f(xk-8/27*h*k1+2*h*k2-3544/2565*h*k3+1859/4104*h*k4-11/40*h*k5,tk+1/2*h,d); 
  x1=xk+25/216*h*k1+1408/2565*h*k3+2197/4104*h*k4-1/5*h*k5;
  x2=xk+16/135*h*k1+6656/12825*h*k3+28561/56430*h*k4-9/50*h*k5+2/55*h*k6;
endfunction

function [t,x,r]=rkf45sev(fsys,x0,ti,tf,tol,hmin,hmax)
//Integración con Runge--Kutta Fehlberg
if ~(isdef('hmin')) then
  hmin=(tf-ti)*1e-9;
end
if ~(isdef('hmax')) then
  hmax=(tf-ti)/50;
end
[f,fc,fe,fi]=fsys();
d=fi(x0,ti);//valor inicial de la variable discreta
fcp=fc(x0,ti,d,-1);
t=zeros(10000);
x=zeros(length(x0),length(t));
x(:,1)=x0;
t(:,1)=ti;
tk=ti;
k=0;
h=hmin; 
r=0;
while tk<tf
  if tk+h>tf then h=tf-tk;end
  k=k+1;
  [x1,x2]=rkf45sev_step(f,x(:,k),tk,d,h); 
  err=norm(x1-x2)/max(norm(x2),%eps)+%eps;

   if (err>tol & h>hmin) then
    //repetimos el paso
    h= (tol/err)^(1/5)*h*0.8;
    k=k-1;
    r=r+1;
   else  
    //antes de aceptar el paso, vemos discontinuidades
      fcn=fc(x2,tk+h,d,-1);

    if or((fcn.*fcp)<0) then
      ind=find(fcn.*fcp<0);
      hMin=0;
      hMax=h;
      dh=h;
      fcpl=fcn;
      //discontinuidad
      while abs(dh)>1e-6
        dh=dh/2;
        h=hMin+dh;
        [x1,x2a]=rkf45sev_step(f,x(:,k),tk,d,h); 
        fcn=fc(x2a,tk+h,d,-1);
        if or(fcn.*fcp<0) then
          ind=find(fcn.*fcp<0);
          hMax=h;
          x2=x2a
          fcpl=fcn;
          //calculamos error por las dudas para el paso siguiente
        else  
          hMin=h;
        end       
      end
      h=hMax;
      fcp=fcpl;
      //detectamos la discontinuidad
      d=fe(x,t,d,ind(1));//ejecutamos el evento
      //aceptamos el paso
      tk=tk+h;
      t(:,k+1)=tk;
      x(:,k+1)=x2;
      h=hmin;
    else 
      fcp=fcn;
      //aceptamos el paso
      tk=tk+h;
      t(:,k+1)=tk;
      x(:,k+1)=x2;
      h= (tol/max(err,0.00000001))^(1/5)*h*0.8;
    end
    
  end
  if h<hmin then h=hmin;end
  if h>hmax then h=hmax;end
end
x=x(:,1:k+1);
t=t(1:k+1);
endfunction
