function [x1,x2]=rkf45_step(f,xk,tk,h)
  k1=f(xk,tk);
  k2=f(xk+1/4*h*k1,tk+1/4*h);
  k3=f(xk+3/32*h*k1+9/32*h*k2,tk+3/8*h); 
  k4=f(xk+1932/2197*h*k1-7200/2197*h*k2+7296/2197*h*k3,tk+12/13*h); 
  k5=f(xk+439/216*h*k1-8*h*k2+3680/513*h*k3-845/4104*h*k4,tk+h); 
  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); 
  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]=rkf45ev(syst,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,nextev]=syst();
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
  //manejo de eventos temporales
  tne=nextev(x(:,k+1),tk)  
  if tk+h>tne then h=tne-tk+1e-12;end
  //  
  k=k+1;
  [x1,x2]=rkf45_step(f,x(:,k),tk,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  
    //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
  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
