function [t,x]=brk34(f,x0,t0,tf,rtol,atol)
   t=zeros(1,10000);
   t(1)=t0; 
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   h=1e-6;
   k=1;
   while t(k)<tf
       if t(k)+h>tf then 
           h=tf-t(k);
       end;
       
       xbrk4=brk4_step(f,x(:,k),t(k),h);
       xrk3=rk3_step(f,x(:,k),t(k),h);
       err=norm(xbrk4-xrk3);
       err0=max(rtol*norm(xrk3),atol);
       if err<err0 then
           x(:,k+1)=xbrk4;
           t(k+1)=t(k)+h;
           k=k+1;
       end
       h=h*(err0/err)^(1/4)*0.8;           
   end
   t=t(1:k);
   x=x(:,1:k);
endfunction


function xk1=rk3_step(f,xk,tk,h)
       k1=f(xk,tk);
       k2=f(xk+h*k1,tk+h);
       k3=f(xk+h*k1/4+h*k2/4,tk+h/2);
       xk1=xk++h*(k1/6+k2/6+2*k3/3);       
endfunction


function xk1=rk4_step(f,xk,tk,h)
       k1=f(xk,tk);
       k2=f(xk+h/2*k1,tk+h/2);
       k3=f(xk+h/2*k2,tk+h/2);
       k4=f(xk+h*k3,tk+h);
       xk1=xk+h*(k1/6+k2/3+k3/3+k4/6);       
endfunction

function xk1=brk4_step(f,xk,tk,h)
       xold=xk;
       J=jacobian(f,xold,tk+h);
       H=-(h*J)^4/24+(h*J)^3/6-(h*J)^2/2+h*J-eye(J);
       xktilde=rk4_step(f,xold,tk+h,-h);
       xnew=xold-pinv(H)*(xk-xktilde);
       while norm(xnew-xold)>1e-8
           xold=xnew;
           J=jacobian(f,xold,tk+h);
           H=-(h*J)^4/24+(h*J)^3/6-(h*J)^2/2+h*J-eye(J);
           xktilde=rk4_step(f,xold,tk+h,-h);
           xnew=xold-pinv(H)*(xk-xktilde);
       end
        xk1=xnew;
endfunction
