function xk1=beulerdae_step(f,xk,tk,h)
  function ff=f_x(z,tk)
    ff=f(z,(z-xk)/h,tk)
  endfunction
  xprev=xk;
  J=jacobiant(f_x,xprev,tk+h);//Jacobiano del Sistema
  invH=pinv(J);//inversa del Hessiano
  xk1=xprev-invH*f_x(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_x,xprev,tk+h);//Jacobiano del Sistema
    invH=pinv(J);//inversa del Hessiano
    xk1=xprev-invH*f_x(xprev,tk+h);//iteracion de Newton
  end  
  
endfunction

function [t,x]=beulerdae(f,x0,ti,tf,h)
//metodo de Backward Euler con iteracion de Newton
//exec('jacobiant.sci')
  t=[ti:h:tf];
  x=zeros(length(x0),length(t));
  x(:,1)=x0;
  for k=1:length(t)-1
   x(:,k+1)=beulerdae_step(f,x(:,k),t(k),h);
  end
endfunction


