function [z,v,info,xxx_newt] = newtonNd(f_newt,x0_newt,Jf_newt,...
                                      E_newt,kmax_newt,dialogo_newt)
//
// Metodo di Newton in R^N.
//
// Approssima uno zero della funzione f_newt: R^N -> R^N con matrice
// jacobiana Jf_newt costruendo la successione x(k).
//
// I valori di f_newt devono essere colonne.
//
// La costruzione della successione si arresta quando:
//   (*) la funzione f_newt si annulla nel punto attuale x(k);
//   (*) la norma del vettore incremento al passo k:
//
//            dx_newt(k) = inv(Jf_newt(x(k))) f(x(k))
//
//       è minore di E_newt (in tal caso, se la successione fosse
//       convergente ad uno zero semplice di f_newt, l'errore assoluto
//       tra z = x(k) e tale zero sarebbe circa pari a E_newt
//       (vedi nota));
//   (*) dopo kmax_newt iterazioni (valore suggerito: kmax_newt = 20).
//
// Se dialogo_newt = "loquace", ad ogni ciclo di iterazione la procedura
// stampa il valore attuale x(k) ed il valore della norma di f_newt(x(k)).
// Disegna infine (nella finestra grafica NFG_newt) due grafici della
// stima dell'errore assoluto.
//
// z: approssimazione finale (zero di f_newt oppure una sua
//    approssimazione);
// v: valore di |f_newt| in z;
// info = 0: individuato valore in cui f_newt si annulla;
//      = 1: criterio di arresto soddisfatto;
//      = 2: la matrice Jf_newt in x(k) è singolare;
//      = 3: superato il numero massimo consentito di iterazioni.
// xxx_newt: matrice contenente la successione di colonne generata x(k).
//
// Se info = 2, si pone z = [] e v = [].
//
info = -1;
NormaScelta = %inf; // Le norme saranno valutate usando questa.
StimaErr = []; // Riga delle stime dell'errore a ciascun passo.
k_newt = 0; // contatore iterazioni effettuate
x_newt = x0_newt;
xxx_newt = x0_newt;
v_f = f_newt(x_newt);
v_J = Jf_newt(x_newt);
if dialogo_newt == "loquace" then
   ListaFinestre_newt = winsid();
   if length(ListaFinestre_newt) ~= 0 then
         FinestraCorrente_newt = gcf(); end;
   NFG_newt = 37; // Numero della Finestra Grafica nella quale la
                  //  procedura disegna.
   clf(37);
   figure_handle_newt = scf(NFG_newt);
   figure_handle_newt.figure_name =...
         "finestra generata da newtonNd.sci (37)";
   printf("x(%d) = [",k_newt);
   for i=1:length(x_newt)-1,
       printf("%6.5e , ",x_newt(i));
   end;
   printf("%6.5e] , |f| = %6.5e\n",x_newt($),norm(v_f,NormaScelta));
end;
//
// ****************** Iterazione
//
while info == -1,
   if v_f == 0 then info = 0; z = x_newt; v = 0; break; end;
   [S,D,P] = lu(v_J); // fattorizzazione EGPP di Jf_newt(x(k))
   if det(D) == 0 then info = 2, z = [], v = []; break; end;
   c = S\(-P*v_f);
   dx_newt = D\c;
   StimaErr($+1) = norm(dx_newt,NormaScelta);
   if StimaErr($) < E_newt then
      info = 1; z = x_newt; v = v_f; break;
   end;
   if k_newt >= kmax_newt then
      info = 3; z = x_newt; v = v_f; break; end;
   xprev_newt = x_newt;
   x_newt = x_newt + dx_newt;
   xxx_newt(:,$+1) = x_newt;
   v_f = f_newt(x_newt);
   v_J = Jf_newt(x_newt);
   k_newt = k_newt+1;
   if dialogo_newt == "loquace" then
      printf("x(%d) = [",k_newt);
      for i=1:length(x_newt)-1,
          printf("%6.5e , ",x_newt(i));
      end;
      printf("%6.5e] , |f| = %6.5e\n",x_newt($),norm(v_f,NormaScelta));
   end;
end;
//
// ****************** Fine iterazione
//
if dialogo_newt == "loquace" then
   select info,
     case 0 then printf("\n\n *** Trovato zero\n\n");
     case 1 then printf("\n\n *** Trovata approssimazione\n\n");
                 printf(" *** Stima errore = %6.5e\n\n", StimaErr($));
     case 2 then
          printf("\n\n *** Errore: Jf(x(%d)) è singolare\n\n",k_newt);
     case 3 then printf("\n\n *** Superato numero max iterazioni\n\n");
   end;
   // Grafici
   StimaErr_x = StimaErr(1:$-1);
   StimaErr_y = StimaErr(2:$);
   subplot(121);
   plot2d("ll",StimaErr_x,StimaErr_y,style=5,frameflag=4);
   xgrid();
   xlabel("stima al passo  k");
   ylabel("stima al passo  k+1");
   title("Stima errore assoluto (scale logaritmiche)","color","blue");
   subplot(122);
   plot(0:length(StimaErr_x)-1,log10(StimaErr_y)./log10(StimaErr_x));
   xgrid();
   title("Stima ordine di convergenza","color","blue");
   xlabel("k");
   ylabel("log10( e(k+1) ) / log10( e(k) )");
   if length(ListaFinestre_newt) ~= 0 then
         scf(FinestraCorrente_newt); end;
end;
endfunction

// Sia f_newt = f di classe C(2) e Z lo zero semplice (ovvero con
// Jf_newt(Z) non singolare) di f che la successione x(k) sta
// approssimando; posto
//
//          e(k) = x(k) - Z   e   F(k) = inv(Jf_newt(x(k))) f(x(k))
//
// si ha:
//
//          e(k) = F(k) + inv(Jf_newt(x(k))) O( |e(k)|^2 )
//
// dunque, quando e(k) è piccolo e la norma di inv(Jf_newt(x(k)))
// non troppo grande, F(k) è una buona stima di e(k).
//
