function [t,x]=BDF34(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;
   for k=1:3
       x(:,k+1)=brk4_step(f,x(:,k),t(k),h);
   end 
   xaux=x(:,1:4);
   k=4;
   while t(k)<tf
       if t(k)+h>tf then 
           h=tf-t(k);
       end;
       xbdf4=BDF4_step(f,xaux,k,t(k),h);
       xrk4=rk4_step(f,x(:,k),t(k),h);
       err=norm(xbdf4-xrk4);
       err0=max(rtol*norm(xbdf4),atol);
       if err<err0 then
           x(:,k+1)=xbdf4;
           xaux(:,k+1)=xbdf4;
           t(k+1)=t(k)+h;
           k=k+1;
           hnew=h*(err0/err)^(1/5)*0.8;
           if hnew>2*h then
               xaux(:,k-3:k)=newval(xaux(:,k-3:k),hnew,h);
               h=hnew;
           end
       else    
           hnew=h*(err0/err)^(1/5)*0.5;
           xaux(:,k-3:k)=newval(xaux(:,k-3:k),hnew,h);
           h=hnew;
       end
   end
   t=t(1:k);
   x=x(:,1:k);
endfunction


function xk1=BDF4_step(f,x,k,tk,h)
       xold=x(:,k);
       s=48/25*x(:,k)-36/25*x(:,k-1)+16/25*x(:,k-2)-3/25*x(:,k-3);
       J=jacobian(f,x(:,k),tk+h);
       H=eye(J)-12/25*h*J;
       xnew=xold-inv(H)*(xold-12/25*h*f(xold,tk+h)-s);
       while norm(xnew-xold)>1e-8
           xold=xnew;
           J=jacobian(f,x(:,k),tk+h);
           H=eye(J)-12/25*h*J;
           xnew=xold-inv(H)*(xold-12/25*h*f(xold,tk+h)-s);
       end
       xk1=xnew;
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


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

M=1/6*[6,0,0,0;11,-18,9,-2;6,-15,12,-3;1,-3,3,-1];
invM=inv(M);
function xnew=newval(xold,hnew,hold)
    dh=hnew/hold;
    H=diag([1,dh,dh^2,dh^3]);
    xnew=(invM*H*M*xold(:,4:-1:1)')';
    xnew=xnew(:,4:-1:1);
endfunction
