function [z,v,info,StimaErr,succ] = newtonNd(f,x0,d1f,E,kmax)
//
// Metodo di Newton per funzioni da R^n in R^n.
//
// Approssima uno zero della funzione f: R^n -> R^n con differenziale
// d1f costruendo una successione x(k) a partire da x0.
//
// La costruzione della successione si arresta quando:
//   (*) la funzione f si annulla nel punto attuale x(k);
//   (*) || x(k+1) - x(k) || risulta minore di E (in tal
//       caso, se la successione fosse convergente ad uno zero semplice
//       di f, l'errore assoluto tra z = x(k) e tale zero sarebbe
//       circa pari ad E);
//   (*) dopo kmax iterazioni (valore predefinito: kmax = 20).
//
//
// z: approssimazione finale (zero di f oppure una sua approssimazione);
// v: valore di f in z;
// info = 0: individuato valore in cui f si annulla;
//      = 1: criterio di arresto soddisfatto;
//      = 2: d1f(x(k)) non invertibile;
//      = 3: superato il numero massimo consentito di iterazioni.
// StimaErr: vettore riga contenente le stime dell'errore calcolate
//           a ciascuna iterazione.
// succ: matrice di colonne x(0), ... ,x(k).
//
//
if ~exists('kmax','l') then kmax = 20; end;
succ = x0;
StimaErr = [];
v_f = f(succ(:,$));
v_d = d1f(succ(:,$));
[S,D,P] = lu(v_d); // fattorizzazione EGPP di d1f(x(k))
if (norm(v_f) ~= 0 & det(D) ~= 0) then passo = -D \ (S \ (P * v_f)); StimaErr = norm(passo); end;
//
// ****************** Iterazione
//
while (norm(v_f) ~= 0 & det(D) ~= 0 & StimaErr($) >= E & (size(succ,'c') - 1) < kmax),
   succ(:,$+1) = succ(:,$) + passo;
   v_f = f(succ(:,$));
   v_d = d1f(succ(:,$));
   [S,D,P] = lu(v_d);
   if det(D) ~= 0 then passo = -D \ (S \ (P * v_f)); StimaErr($+1) = norm(passo); end;
end;
//
// ****************** Fine iterazione
//
z = succ(:,$); v = v_f; info = 1;
if norm(v_f) == 0 then info = 0;
   else if det(D) == 0 then info = 2;
           else if (size(succ,'c') - 1) == kmax then info = 3; end;
        end;
end;
endfunction
//
