function [t,x]=rk34(f,x0,ti,tf,rtol)
  n=length(x0);
  t=zeros(1,10000);
  x=zeros(n,10000);
  x(:,1)=x0;
  t(1)=ti;
  k=1;
  h=1e-8;
  hmax=(tf-ti)/20;
  hmin=1e-12;
  atol=1e-9;
  while t(k)<tf
    if t(k)+h>tf then 
      h=tf-t(k); 
    end
    xa=rk3_step(f,x(:,k),t(k),h);
    xb=rk4_step(f,x(:,k),t(k),h);
    e=norm(xa-xb);
    e0=max(rtol*norm(xb),atol);
    if e<=e0 then
      //aceptamos el paso
      x(:,k+1)=xb; //tiene mejor estabilidad
      t(k+1)=t(k)+h;
      k=k+1;
    end
    h=h*(e0/e)^(1/4)*0.8;
    h=min(h,hmax);
    h=max(h,hmin);
  end
  t=t(1:k);
  x=x(:,1:k);
endfunction
