function [t,x]=BDF4(f,x0,t0,tf,h)
   t=[t0:h:tf];
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   
   for k=1:3
       x(:,k+1)=brk4_step(f,x(:,k),t(k),h);
   end 
   
   for k=4:length(t)-1
       x(:,k+1)=BDF4_step(f,x,k,t(k),h);
   end 
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
