// determiniamo l'errore massimo nella soluzione
// dell'equazione differenziale y'=0.001y(1000-y)
// con l'algoritmo di Eulero
// nell'intervallo 0<t<20
// Impiega alcuni minuti....
float f(float t, float y) {
   float f=0.001*y*(1000-y);
   return f;
}
float soluzione(float t) {
   float s=1000./(exp(-t)*999+1);
   return s;
}
void RKErrore1() {
   // TTree che conterra il risultato delle simulazione
   // Apro file
    TFile * file=new TFile("RKErr1.root","RECREATE");
   // Creo il TTree
    TTree * T= new TTree("T","Equazioni Differenziali");
   // Le variabili del TTree
    float t=0,y=0,y2=0,y4=0,s=0;
    float NPoint,step,errEU,errRK2,errRK4;
    T->Branch("NPoint",&NPoint,"NPoint/F");
    T->Branch("step",&step,"step/F");
    T->Branch("errEU",&errEU,"errEU/F");
    T->Branch("errRK2",&errRK2,"errRK2/F");
    T->Branch("errRK4",&errRK4,"errRK4/F");
//step usato per la discretizzazione
    float h;
    float localerr, localerr2, localerr4;
// variamo opportunamente lo step
// aumentando il numero di punti in cui
// calcolo la soluzione
    int N=20;
    NPoint = 10.;
    for (int n=0; n<N;n++) {
        NPoint=10*pow(2,n);
        h=20./NPoint;
        step=h;
        errEU=-10;
        errRK2=-10;
        errRK4=-10;
// condizioni iniziali
        t=0;
	y =1;  /// soluzione metodo di Eulero
        y2=1;  /// soluzione metodo di RK del secondo ordine 
        y4=1;  /// soluzione metodo di RK del  quarto ordine
// Loop principale
        float K1,K2,K3,K4;
	localerr = 0.;
	localerr2 = 0.;
	localerr4 = 0.;
        for (int i=0; i<NPoint; i++)  {

// formula di Runge Kutta del primo  ordine = Eulero 
	    K1=f(t,y);
            y=y+h*K1;
	  
// formula di Runge Kutta del secondo ordine
            K1=f(t    , y2       );
            K2=f(t+h/2, y2+h/2*K1);
            y2=y2+h*K2;

// formula di Runge Kutta del quarto ordine
            K1 = f(t    , y4        );
            K2 = f(t+h/2, y4 + h * K1/2.);
            K3 = f(t+h/2, y4 + h * K2/2.);
            K4 = f(t+h  , y4 + h * K3   );
            y4 = y4 + h * ( K1 + 2*K2 + 2*K3 + K4 ) / 6.;

	    // aggiorno t 
            t=t+h;
	  
// Valore esatto
            s=soluzione(t);
// Calcolo dell'errore massimo
            localerr=s-y;
            if (fabs(localerr) >errEU )errEU =fabs(localerr);
            localerr2=s-y2;
            if (fabs(localerr2)>errRK2)errRK2=fabs(localerr2);
            localerr4=s-y4;
            if (fabs(localerr4)>errRK4)errRK4=fabs(localerr4);

        }
	std::cout<<"Npoint = "<< NPoint<<" step = " <<h<<" errMax = "<<errEU<<" "<<errRK2<<" "<<errRK4<<std::endl;
        T->Fill();

    }
//
    T->SetMarkerStyle(20);
    T->Write();
    file->Close();

}