function [x,M]=broyden_g(func,x,M)

% Solving f(x)=0 with the Broyden's (good) method

tol=10^(-10); k=1; kmax=500;
f=feval(func,x); err=norm(f); dx=-M*f; d=1;

while err>tol ,
   xn=x+dx; fn=feval(func,xn); errn=norm(fn);
   z=M*fn; s=1/(dx'*(z+dx)); w=s*M'*dx;
   M=M-z*w';
   if errn<err , x=xn; dx=-(1-w'*fn)*z; err=errn; f=fn;
   else , dx=-M*f; end; % We don't take the bad step
   k=k+1; if k>kmax , err, err=0; x=NaN*x; end; end


   
