function [x,nit,flag]=broydeninv(F,H,x)
nmax=1000;
toll=1.e-12;
flag=1;
Fx=feval(F,x);
for nit=1:nmax
    s=-H*Fx;
    x=x+s;
    inc=norm(s); 
    if abs(inc)<toll, flag=0; return, end
    Fx=feval(F,x);
    z=H*Fx;
    beta=s'*s+s'*z;
    H=H-1/beta*z*s'*H;
end