function [T,X,PASSO] = RK12_pv(x0,t0,tf,F,E,LAMBDA,HMIN,TAU)
//
// 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 esplicito (RK di ordine 1) - a passo
// variabile - affiancato, per la scelta del passo, dal metodo di
// Heun (RK di ordine 2).
//
// x0: condizione iniziale (colonna di n elementi)
// t0: istante iniziale (numero reale)
// tf: istante finale (numero reale)
// F: function che definisce l'equazione differenziale - F(t,x) deve
//    essere una colonna di n numeri reali
// E: valore massimo della stima dell'errore locale (numero reale)
// LAMBDA: numero reale che stabilisce il valore massimo del passo
//         (OPZIONALE - valore predefinito: 1d-5)
// HMIN: valore minimo consentito del passo
//       (OPZIONALE - valore predefinito: (tf - t0) / 1d6)
// TAU: valore del passo per il calcolo delle stime utilizzate
//      nel calcolo di h(k) (OPZIONALE - valore predefinito: (tf - t0) / 1d3)
//
// T = [t(0),...,t(N)]: riga contenente gli istanti di integrazione
// X = [x(0),...,x(N)]: matrice n x (N+1) contenente le approssimazioni
// PASSO = [h(0),...,h(N-1)]: riga contenente i passi di integrazione
//
// Valore degli argomenti opzionali
//
if ~exists('LAMBDA','l') then LAMBDA = 1d-5; end;
if ~exists('HMIN','l') then HMIN = (tf - t0) / 1d6; end;
if ~exists('TAU','l') then TAU = (tf - t0) / 1d3; end;
//
// Inizializzazione delle variabili di uscita
//
T(1,1) = t0;
X(:,1) = x0;
PASSO = [];
//
// ciclo principale
//
while (T(1,$) < tf), // arresta la costruzione se ha raggiunto tf
  //
  // scelta del passo
  //
  // XX1 = X(:,$) + F(T(1,$),X(:,$)) * TAU;
  ST1 = F(T(1,$),X(:,$));
  ST2 = F( T(1,$) + TAU,X(:,$) + ST1 * TAU );
  // XX2 = X(:,$) + ( (ST1 + ST2)/2 ) * TAU;
  //
  //     XX1 - XX2 = (ST1 - ST2)/2 * TAU
  //
  d = max(LAMBDA, norm( ((ST1 - ST2)/2) * TAU ));
  PASSO(1,$+1) = min(sqrt(E/d) * TAU, tf - T(1,$));
  //
  // calcolo approssimazione e nuovo istante di integrazione
  //
  X(:,$+1) = X(:,$) + F(T(1,$),X(:,$)) * PASSO(1,$);
  T(1,$+1) = T(1,$) + PASSO(1,$);
  //
  // arresta la costruzione se il passo calcolato risulta troppo
  // piccolo e non ha raggiunto tf
  //
  if (PASSO(1,$) < HMIN) & (T(1,$) < tf) then break; end;
  //
end;
//
// Verifica se l'integrazione ha raggiunto tf
//
if T(1,$) < tf then
   printf("\n\nIntegrazione interrotta a T = %3.2e", T(1,$));
end;
//
endfunction
//
// Esempio per assegnare valori ai parametri opzionali:
//
//    [T,X,PASSO] = TS_1_pv(x0,t0,tf,F,G,E,HMIN = y);
//
//         => LAMBDA = valore predefinito, HMIN = y, TAU = valore predefinito
//
