// predict_correct.cpp : Defines the entry point for the console application.
//


/*
*   ADAMS-FORTH ORDER PREDICTOR-CORRECTOR ALGORITHM 5.4
*
*   To approximate the solution of the initial value problem
*          y' = f(t,y), a <= t <= b, y(a) = alpha,
*   at N+1 equally spaced points in the interval [a,b].
*
*   INPUT:   endpoints a,b; initial condition alpha; integer N.
*
*   OUTPUT:  approximation w to y at the (N+1) values of t.
*/

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

#define true 1
#define false 0
#define MAX_N 2000

double F(double, double);
void RK4(double , double , double , double *, double *);

int main( )
{
      double T[MAX_N],W[MAX_N], Wp;
      double A,B,ALPHA,H,T0,W0;
      int I,N,J;

      /* STEP 1 */
	  B = 2.0;
	  A= 0.0;
	  ALPHA = 0.5;
	  N = 10;
      H = (B - A) / N;
      T[0] = A;
      W[0] = ALPHA;
      printf("T[%d] = %5.3f; W[%d]= %11.7f\n", 0, T[0], 0, W[0]);
      
      /* compute starting values using Runge-Kutta method */
      /* given in a subroutine */
	  for (I=1; I<=3; I++)
	  {
         RK4(H, T[I-1], W[I-1], &T[I], &W[I]);
         printf("T[%d] = %5.3f; W[%d] = %11.7f\n", I, T[I], I, W[I]);
      }

	  /* Prediction-correction step*/
      for (I=4; I<=N; I++) 
	  {
         /* T0, W0 will be used in place of t, w resp. */
         T[I] = A + I * H;
         /* predict W(I) */
         Wp = W[I-1]+H*(55.0*F(T[I-1],W[I-1])-59.0*F(T[I-2],W[I-2])
              +37.0*F(T[I-3],W[I-3])-9.0*F(T[I-4],W[I-4]))/24.0;

         /* correct W(I) */
         W[I] = W[I-1]+H*(9.0*F(T[I],Wp)+19.0*F(T[I-1],W[I-1])
              -5.0*F(T[I-2],W[I-2])+F(T[I-3],W[I-3]))/24.0;
         
         printf("T[%d] = %5.3f; W[%d] = %11.7f\n", I, T[I], I, W[I]);
         
      }

	  return 0;
}


/*  Change function F for a new problem    */
double F(double T, double Y)
{
   double f; 

   f = Y - T*T + 1.0;
   return f;
}

void RK4(double H, double T0, double W0, double *T1, double *W1)
{
   double K1,K2,K3,K4;
   
   *T1 = T0 + H;
   K1 = H * F(T0, W0);
   K2 = H * F(T0 + 0.5 * H, W0 + 0.5 * K1);
   K3 = H * F(T0 + 0.5 * H, W0 + 0.5 * K2);
   K4 = H * F(*T1, W0 + K3);
   *W1 = W0 + (K1 + 2.0 * (K2 + K3) + K4) / 6.0;
}