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