function xk1=bdf3dae_step(f,xk,xkm1,xkm2,tk,h)
  function ff=f_xxx(z,tk)
    ff=f(z,(11/6*z-3*xk+3/2*xkm1-1/3*xkm2)/h,tk)
  endfunction
  xprev=xk;
  J=jacobiant(f_xxx,xprev,tk+h);//Jacobiano del Sistema
  invH=pinv(J);//inversa del Hessiano
  xk1=xprev-invH*f_xxx(xprev,tk+h);//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_xxx,xprev,tk+h);//Jacobiano del Sistema
    invH=pinv(J);//inversa del Hessiano
    xk1=xprev-invH*f_xxx(xprev,tk+h);//iteracion de Newton
  end  
  
endfunction

function [t,x]=bdf3dae(f,x0,ti,tf,h)
//metodo de BDF3 con iteracion de Newton
  t=[ti:h:tf];
  x=zeros(length(x0),length(t));
  [tin,xin]=radau3(f,x0,ti,ti+2*h,h);
  x(:,1:3)=xin;
  for k=3:length(t)-1
   x(:,k+1)=bdf3dae_step(f,x(:,k),x(:,k-1),x(:,k-2),t(k),h);
  end
endfunction


