function [t,x]=brk34(f,x0,t0,tf,reltol)
    abstol=1e-9;
    hmin=1e-9*(tf-t0);
    hmax=(tf-t0)/50;   
    h=hmin;
    N=3;
    t=zeros(1,100000);
    x=zeros(length(x0),length(t));
    t(1)=t0;
    x(:,1)=x0;
    k=1;
    while t(k)<tf
        x1=rk3_step(f,x(:,k),t(k),h);
        x2=brk4_step(f,x(:,k),t(k),h);
        e=norm(x1-x2)+1e-12;
        tol=reltol*norm(x1)+abstol;
        if e<tol then
            x(:,k+1)=x2;
            t(k+1)=t(k)+h;
            k=k+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]=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);
    x1=xk+h*(k1/6+k2/6+k3*2/3);
endfunction

function xk1=brk4_step(f,xk,tk,h)
   xk1=xk;//initial guess
   xk_l=rk4_step(f,xk1,tk+h,-h);
   it=0;
   while norm(xk_l-xk)>1e-6*(norm(xk)+1e-3)&it<10
       J=jacobiant(f,xk1,tk+h);
       invH=pinv(eye(J)-J*h+(J*h)^2/2-(J*h)^3/6+(J*h)^4/24);
       xk1=xk1-invH*(xk_l-xk)
       xk_l=rk4_step(f,xk1,tk+h,-h);
       it=it+1;
    end
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/6*(k1+2*k2+2*k3+k4);
endfunction

function J=jacobiant(f,x,t)
 //estimates the Jacobian df/dx evaluated at (x,t) 
  n=length(x);
  J=zeros(n,n);
  fx=f(x,t);
  for i=1:n
    x1=x;
    dx=0.00001*x1(i)+1e-10;
    x1(i)=x1(i)+dx;
    J(:,i)=(f(x1,t)-fx)/dx;
  end
endfunction  
            
    