function [t,x]=beuler_dae(f,x0,t0,tf,h)
   t=[t0:h:tf];
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   for k=1:length(t)-1
       x(:,k+1)=beuler_dae_step(f,x(:,k),t(k),h);
   end 
endfunction

function xk1=beuler_dae_step(f1,xk,tk,h)

       function res=ftilde(xk1,tk1)
         res=f1(xk1,(xk1-xk)/h,tk1)
       endfunction
       
       xold=xk;
       H=jacobian(ftilde,xold,tk+h);
       xnew=xold-inv(H)*ftilde(xold,tk+h);
       while norm(xnew-xold)>1e-8
           xold=xnew;
           H=jacobian(ftilde,xold,tk+h);
           xnew=xold-inv(H)*ftilde(xold,tk+h);
       end
       xk1=xnew;
endfunction
