/* Adams Fourth Order Predictor-Corrector Program 
 * This solves initial value problems using the Adams Predictor-Corrector
 * Algorithm. It inputs The range of interval of time steps and the initial
 * value alpha where y(0) = alpha.
 *
 * apc V. 0.9 written in C by Stefan Joe-Yen,  Nov 1995
 */

/* Standard Includes */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Define the function being solved y' = f(t,y)  and the number of time 
steps N */ 

#define f(t,y) ((float) (y - (t*t) + 1))
#define N 10

main()
{
   float a, b, h, t[N+1], w[N+1], alpha, K1, K2, K3, K4;
   int i,j;

   /* Input the endpoints and initial value */
   printf("please enter a, b, and alpha>");
   scanf("%f %f %f", &a, &b, &alpha);

   /* Set up Initial Conditions */
   h = (b - a)/N;
   t[0] = a;
   w[0] = alpha;
   printf("t[0] = %f, w[0] = %f\r\n", t[0], w[0]);

   /* For the first 3 Time steps use Runge-Kutta to estimate solutions */

   for (i=1; i<=3; i++) {
      K1 = h * f(t[i-1],w[i-1]);
      K2 = h * f(t[i-1] + h/2,w[i-1] + K1/2);
      K3 = h * f(t[i-1] + h/2,w[i-1] + K2/2);
      K4 = h * f(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("t[%d] = %f, w[%d] = %f\r\n", i, t[i], i, w[i]);
   }

   /* For subsequent time steps predict and correct the estimated 
      solution */

   for (i=4; i<=N; i++) {
      t[i] = a + i * h;

      /* predict w[i] */
      w[i] = w[3] + h*((55 * f(t[3],w[3])) - (59 * f(t[2], w[2])) + (37 * 
                       f(t[1],w[1])) - (9 * f(t[0], w[0])))/24;

      /* Correct w[i] */
      w[i] = w[3] + h*((9 * f(t[i],w[i])) + (19 * f(t[3],w[3])) - (5 * 
                       f(t[2],w[2])) + f(t[1],w[1]))/24;

      /* Output estimated solution at time step t */
      printf("t[%d] = %f, w[%d] = %f\r\n", i, t[i], i, w[i]);
      
     
     /* Shift values to prepare for the next iteration */
      for (j=0; j<=2; j++) {
         t[j] = t[j+1];
         w[j] = w[j+1];
      }
      t[3] = t[i];
      w[3] = w[i];
   }

}