function [t,x]=beuler(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_step(f,x(:,k),t(k),h);
   end 
endfunction

function xk1=beuler_step(f,xk,tk,h)
       xold=xk;
       J=jacobian(f,xk,tk+h);
       H=h*J-eye(J);
       xnew=xold-inv(H)*(xk+h*f(xold,tk+h)-xold);
       while norm(xnew-xold)>1e-8
           xold=xnew;
           J=jacobian(f,xold,tk+h);
           H=h*J-eye(J);
           xnew=xold-inv(H)*(xk+h*f(xold,tk+h)-xold);
       end
       xk1=xnew;
endfunction
