function xk1=brk4_step(f,xk,tk,h)
  //armamos la iteracion de Newton
  I=eye(length(xk),length(xk))//matriz identidad
  xprev=xk;
  J=jacobiant(f,xprev,tk+h);//Jacobiano del Sistema
  invH=pinv(I-h*J+(h*J)^2/2-(h*J)^3/6+(h*J)^4/24);//inversa del Hessiano
  F=rk4_step(f,xprev,tk+h,-h)-xk;//función que debe ser cero
  xk1=xprev-invH*F;//iteracion de Newton
  l=0;
  while norm(xk1-xprev)>1e-6*(norm(xprev)+1e-3)&l<20
    l=l+1;
    xprev=xk1;
    J=jacobiant(f,xprev,tk+h);//Jacobiano del Sistema
    invH=pinv(I-h*J+(h*J)^2/2-(h*J)^3/6+(h*J)^4/24);//inversa del Hessiano
    F=rk4_step(f,xprev,tk+h,-h)-xk;//función que debe ser cero
    xk1=xprev-invH*F;//iteracion de Newton
  end    
endfunction

function xk1=bdf4_step(f,xk,xkm1,xkm2,xkm3,tk,h)
  //armamos la iteracion de Newton
  I=eye(length(xk),length(xk))//matriz identidad
  xprev=xk;
  J=jacobiant(f,xprev,tk+h);//Jacobiano del Sistema
  invH=pinv(12/25*h*J-I);//inversa del Hessiano
  F=12/25*h*f(xprev,tk+h)-xprev+48/25*xk-36/25*xkm1+16/25*xkm2-3/25*xkm3;//función que debe ser cero
  xk1=xprev-invH*F;//iteracion de Newton
  l=0;
  while norm(xk1-xprev)>1e-6*(norm(xprev)+1e-3)&l<20
    l=l+1;
    xprev=xk1;
    J=jacobiant(f,xprev,tk+h);//Jacobiano del Sistema
    invH=pinv(12/25*h*J-I);//inversa del Hessiano
    F=12/25*h*f(xprev,tk+h)-xprev+48/25*xk-36/25*xkm1+16/25*xkm2-3/25*xkm3;//función que debe ser cero
    xk1=xprev-invH*F;//iteracion de Newton
  end    
endfunction

function [t,x]=bdf4(f,x0,ti,tf,h)
//metodo de BDF4 con iteracion de Newton
//exec('jacobiant.sci')
//exec('rk4.sci')
t=[ti: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),x(:,k-1),x(:,k-2),x(:,k-3),t(k),h);
end


endfunction


