function [T,X,PASSO] = eulero_pv(x0,t0,tf,fct,EL_MAX,dialogo)
//
//  Integra numericamente, sull'intervallo [t0,tf], il problema
//  di Cauchy in R(n):
//  .
//  x = F(t,x)
//  x(t0) = x0
//
//  con il metodo di Eulero in avanti a passo variabile.
//
//  x0: condizione iniziale (colonna di n elementi)
//  t0: istante iniziale
//  tf: istante finale
//  fct: function per F - fct(t,x) deve essere una colonna
//  EL_MAX: errore locale massimo consentito
//  dialogo: se "loquace" mostra gli istanti di integrazione e, ad ogni iterazione,
//           se è stato necessario aumentare o diminuire la lunghezza del passo
//
//  T = [T(1),...,T(N)], nodi
//  X: matrice n x N - la colonna X(:,i) è la soluzione numerica
//     all'istante T(i)
//  PASSO: riga con PASSO(k) = h tale che T(k+1) = T(k) + h
//
n = length(x0); // determina il numero di equazioni del sistema
h_min = (tf - t0)/1d7; // passo minimo consentito
h_max = (tf - t0)/10;  // passo massimo consentito
T = [];
X = [];
PASSO = [];
//
T(1,1) = t0;
X(:,1) = x0;
//
h = h_max/10;  // passo iniziale
//
// ciclo principale
//
while (T(1,$) < tf) & (h >= h_min),  // l'iterazione si arresta se
                                     // si è raggiunto tf o se non si
                                     // è riuscito a rendere la stima dell'errore
                                     // locale inferiore a EL_MAX
    h_max_loc = min(tf - T(1,$), h_max);
    if h > h_max_loc then h = h_max_loc; end;
    xx1 = X(:,$) + h*fct(T(1,$),X(:,$));
    // stima EL
    xx2 = X(:,$) + (h/2)*fct(T(1,$),X(:,$));
    xx2 = xx2 + (h/2)*fct(T(1,$)+h/2,xx2);
    StimaEL = 2*norm(xx1 - xx2);
    // fine stima EL
    // decide se accettare il passo in base al valore di StimaEL
    if StimaEL > EL_MAX then  // passo troppo lungo: riduci h...
       h = h/2; if dialogo == "loquace" then printf("-"); end;
       elseif (StimaEL < EL_MAX/4) & (h < h_max_loc) then
              // passo troppo corto: aumenta h...
              h = min(1.18*h, h_max_loc); // 1.18 ~= 2^(1/4): 4 cicli per
                                          // "recuperare il dimezzamento"
              if dialogo == "loquace" then printf("+"); end;
              else T(1,$+1) = T(1,$) + h;  // valore di h accettato...
                   if dialogo == "loquace" then printf("\nT = %3.2e",T($)); end;
                   X(:,$+1) = xx1;
                   PASSO(1,$+1) = h;
    end;
end;
if dialogo == "loquace" then printf("\n"); end;
//
