function [x,res,niter] = newtonsys(Ffun,Jfun,x0,tol,...
                                nmax, varargin)
%NEWTONSYS calcola una radice di un sistema non lineare
%   [ZERO,RES,NITER]=NEWTONSYS(FFUN,JFUN,X0,TOL,NMAX)
%   calcola il vettore ZERO, radice di un sistema non
%   lineare definito nella function FFUN con matrice
%   Jacobiana definita nella function JFUN a partire
%   dal vettore X0. RES contiene il valore del residuo
%   in ZERO e NITER il numero di iterazioni necessarie
%   per calcolare ZERO. FFUN e JFUN sono function
%   definite tramite M-file

niter = 0;
err = tol + 1;
x = x0;
while err >= tol & niter < nmax
    J = Jfun(x,varargin{:});
    F = Ffun(x,varargin{:});
    delta = - J\F;
    x = x + delta;
    err = norm(delta);
    niter = niter + 1;
end
res = norm(Ffun(x,varargin{:}));
if (niter==nmax & err> tol)
    fprintf(['Il metodo non converge nel massimo ',...
       'numero di iterazioni. L''ultima iterata\n',...
       'calcolata ha residuo relativo pari a %e\n'],F);
else
    fprintf(['Il metodo converge in %i iterazioni',...
            ' con un residuo pari a %e\n'],niter,F);
end
return
