/* AMB.c   -- solution of IVP by the the Adams-Bashforth-Moulton
 Predictor-Corrector scheme.
 this program provides an approximate (numerical) solution to
 the well-posed initial value problem y' = f(t,y), a <= t <= b, (1)
 y(a) = alpha at (n+1) equally-spaced abscissa points in the interval
 [a,b].  note that y is taken as the dependent function of the
 independent variable t.

 The ABM Predictor-Corrector Method constructs an approximation w[i]
 to the true solution
 y(t[i]) at each mesh point in the interval [a,b] for each
 i = 1,2, ... , n by the recursive scheme

 w[0] = alpha
 w[i+1] = w[i] + (K1 + 2K2 + 2K3 + K4)/6
 The above is the same as the RK4 solution which is used as the
 Predictor.  This solution is used for the first 3 node points.
 The succeeding node points are calculuated using the "corrector"
 scheme.  */

#include <stdio.h>
#include <math.h>

double func(double t,double y);      //function prototype for y ' = f(t,y)

int main(void)
{
   double a,b,alpha,h,t[4],tcurr,w[4],wpred,wcorr,K1,K2,K3,K4;
   int n,i,j;
   printf("Adams-Bashforth-Moulton Predictor-Corrector\nMethod for solution of y' = f(t,y) on [a,b], y(a)=alpha\n");
   printf("Please enter interval limits [a,b], b > a:\n");
   scanf("%lf %lf",&a,&b);
   while (b<=a)
   {
      printf("b must be greater than a.  Please re-enter a and b: \n");
      scanf("%lf %lf",&a,&b);
   }
   printf("Please enter number of steps n to take in [a,b]:\n");
   scanf("%d",&n);
   while (n<=0)
   {
	  printf("n must be a positive integer. Please re-enter n: \n");
      scanf("%d",&n);
   }

   h = (b-a)/n;
   printf("Please enter initial condition alpha, y[%lf] = alpha:\n",a);
   scanf("%lf",&alpha);
   /*  INITIALIZATION       */
   i = 0;
   t[0]=a;
   w[0] = alpha;
   /* END OF INITIALIZATION */
   printf("\n\nIteration\t\tt\t\t\ty\n");
   printf("%d\t\t\t%6.6lf\t\t%6.6lf\n",i,t,w);
   for (i = 1;i<=3;i++)
   {
	   K1 = h*func(t[i-1],w[i-1]);
	   K2 = h*func(t[i-1]+h/2,w[i-1]+K1/2);
	   K3 = h*func(t[i-1]+h/2,w[i-1]+K2/2);
	   K4 = h*func(t[i-1]+h,w[i-1]+K3);
	   w[i] = w[i-1] + (K1 + 2*K2 + 2*K3 + K4)/6.;
	   t[i] = a + i*h;
	   printf("%d\t\t\t%6.6lf\t\t%6.6lf\n",i,t[i],w[i]);
   }

   for (i = 4;i<=n;i++)
   {
	  tcurr = a+i*h;
	  wpred = w[3]+h*(55*func(t[3],w[3])-59*func(t[2],w[2])+37*func(t[1],w[1])-9*func(t[0],w[0]))/24; //predict w
	  wcorr = w[3]+h*(9*func(tcurr,wpred)+19*func(t[3],w[3])-5*func(t[2],w[2])+func(t[1],w[1]))/24; //correct w
	  printf("%d\t\t\t%6.6lf\t\t%6.6lf\n",i,tcurr,wcorr);
      for(j=0;j<=2;j++)
      {
		  t[j] = t[j+1];
		  w[j] = w[j+1];
	  }
	  t[3] = tcurr;
	  w[3] = wcorr;

   }

   return 0;
}

double func(double t,double y)
{
	return -y + t + 1;
}