//************************************************************************ 
//  Solve the IVP   
//  x'=t^2-x, x(0)=1
//  over 0<t<5 with the Adams-Bashforth-Moulton predictor-corrector method.
//  The three starting values x1, x2, x3 can be estimated using the classical RK4
//  method.
//  The exact solution is
//
//  x(t)=t^2-2*t+2-exp(-t).
//************************************************************************
/*
compile with	
gcc -o abm abm.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/


#include "/usr/local/dislin/dislin.h"
#include <math.h>
#include <stdio.h>
#include<stdlib.h>



//define functions for r.h.s.
double fun(double t, double x) {
  return pow(t,2)-x;
}
//
int step=0.0;
int i=0;
int k=3;
double max_steps = 5.0; // final time
//double N=2000.0; // amount of steps

//
float T[2000];
float X[2000];
float xe[2000];

double t = 0.0; //initial conditions
double x = 1.0;
float p;
double h;
double h1;


//
//
double k1, k2, k3, k4;

double f0, f1, f2, f3, f4; /* Function values                */

main ()
{ 
T[0]=0.;
X[0]=1.;

h=0.05;
// Calculate starting values x_i, i=1,2,3 with the classical RK4 method
 for (i = 0; i<=3; i++)
  {
   t  = T[i];
   x  = X[i];
   k1 = h * fun(t, x);
    
    
    k2 = h * fun(t + 0.5 * h, x + 0.5 * k1);
    
    k3 = h * fun(t + 0.5 * h, x + 0.5 * k2);
    
    k4 = h * fun(t + h, x+ k3);
    
    X[i+1]= x+(k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
    
    T[i+1] = t + h;
  }
    f0 = fun(T[0],X[0]);  
    f1 = fun(T[1],X[1]);
    f2 = fun(T[2],X[2]);
    f3 = fun(T[3],X[3]);
    h1=h/24.0;

 while (T[k] < max_steps)
   {    
        /* Predictor */
        p = X[k] + h1 * (-9.0*f0 + 37.0*f1 - 59.0*f2 + 55.0*f3);
        /* Next time point */
        T[k+1] =h*(k+1);
        /* Evaluate f(t,y) */
        f4 = fun(T[k+1],p);
        /* Corrector */
        X[k+1] = X[k] + h1 * (f1 - 5.0*f2 + 19.0*f3 + 9.0*f4);
        /* Update the values */
        f0 = f1;
        f1 = f2;
        f2 = f3;
        f3 = fun(T[k+1], X[k+1]);
        k++;
   }
  // show calculated values
  for ( i = 0; i <= k; i++)
    {
       printf("i = %d, T[i] = %lf, X[i] = %lf\n", i, T[i], X[i]);
    }

// plot resulting solution
  int nn=i-1;
  winsiz (600, 600);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini();
  pagera();
  texmod ("on");
  axspos(450,1800);
  axslen(1800,1200);

  name("t","x");
  name("x(t)","y");

  graf(0.,max_steps,0.,1.0,0.0,18.0,0.,1.0);
  height(50);
  title();

  color("red");
  curve(T,X,nn);
// plot exact solution
for (i=0; i<nn; i++)
{
xe[i]=pow(T[i],2)-2*T[i]+2-exp(-T[i]);
}

color("blue");
incmrk(-5);
curve(T,xe,nn);

  endgrf();
  disfin();
}