/*
*
*   To approximate the solution of the initial value problem
*          y' = 0, 0 <= t <= b, y(0) = 1,
*   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 5000

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 = 10.0;
	  A= 0.0;
	  ALPHA = 9.4;
	  N = 30;
      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]);
      
       // RK method to get W1.
      for(I = 1; I <=1; I++)
      {
          RK4(H, T[0], W[0], &T[1], &W[1]);
          printf("T[%d] = %5.3f; W[%d]= %11.7f\n", I, T[I], I, W[I]);
      } 

      /* An unstable 2-step  method*/
      for (I=2; I<=N; I++) 
      {
         /* T0, W0 will be used in place of t, w resp. */
         T[I] = A + I * H;

         /* correct W(I) */
         W[I] = -4.0*W[I-1]+ 5.0*W[I-2] + H*(4.0*F(T[I-1], W[I-1])+2.0*F(T[I-2],W[I-2]));
         
         printf("T[%d] = %5.3f; W[%d] = %14.13f\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 = 0.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;
}