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

/*
*    Adams-Bashforth four-step explicit Multistep method.
*    RUNGE-KUTTA (ORDER 4) ALGORITHM is used to compute w1
*
*    TO APPROXIMATE THE SOLUTION TO THE INITIAL VALUE PROBLEM:
*               Y' = F(T,Y), A<=T<=B, Y(A) = ALPHA,
*    AT (N+1) EQUALLY SPACED NUMBERS 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 "stdafx.h"

#include<stdio.h>
#include<math.h>
#include<iostream>
#include <fstream>
#include <string>

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

static   double F(double, double, double);
static   double F2(double, double, double);

using namespace std;

int main()
{
      string str[30][3];
      ofstream myoutput("output.txt");

      double A,B,ALPHA,H,T,W[MAX_N][2], K1,K2,K3,K4,L1,L2,L3,L4;
      int I,N;

      A = 0.0;
      B = 1.0;
      N = 100;

      cout.setf(ios::fixed,ios::floatfield);
      cout.precision(5);

      /* STEP 1 */
      H = (B - A) / N;
      T = A;
      W[0][0] = 5.0; // initial value
      W[0][1] = -10.0;

      /* STEP 2 --- Use order 4 RK method to get w1, w2, w3 */
      /// NOTE: The "for" loop starts with I = 1.
      if (myoutput.is_open())
    {
      //For F1 and F2
      for (I=1; I<=2; I++)
      {
         /* STEP 3 */
         /* Compute K1, K2 RESP. */
         K1 = H*F(T, W[I-1][0],W[I-1][1]);
         L1 = H*F2(T, W[I-1][0],W[I-1][1]);
         K2 = H*F(T + H/2.0, W[I-1][0] + K1/2.0, W[I-1][1] + L1/2.0);
         L2 = H*F2(T + H/2.0, W[I-1][0] + K1/2.0, W[I-1][1] + L1/2.0);
         //K3 = H*F(T + H/2.0, W[I-1] + K2/2.0);
         //K4 = H*F(T + H, W[I-1] + K3);
         K3 = H*F(T + H, W[I-1][0] + 2.0*K2 - K1, W[I-1][1] + 2.0*L2 - L1);
         L3 = H*F2(T + H, W[I-1][0] + 2.0*K2 - K1, W[I-1][1] + 2.0*L2 - L1);

         /* STEP 4 */
         /* COMPUTE W(I) */
         //W[I] = W[I-1] + 1/6.0*(K1 + 2.0*K2 + 2.0*K3 + K4);
         W[I][0] = W[I-1][0] + 1/6.0*(K1 + 4.0*K2 +K3);
         W[I][1] = W[I-1][1] + 1/6.0*(L1 + 4.0*L2 +L3);


         /* COMPUTE T(I) */
         T = A + I * H;

         /* STEP 5 */
         cout <<"F1 - F2: At time "<< T <<" the solution = "<< W[I][0] << "," << W[I][1] << endl;

             for (int b=0;b<2;b++)
             {
                 str[I-1][b] = W[I][b];
                 myoutput << str[I-1][b]<< " ";
             }
      }

      /* STEP 6---Use Adam-Bashforth 4-step explicit method */
      //for(I = 4; I <= N; I++)
      for(I = 3; I <= N; I++)
      {
          //K1 = 55.0*F(T, W[I-1]) - 59.0*F(T-H, W[I-2]) + 37.0*F(T-2.0*H, W[I-3]) - 9.0*F(T-3.0*H, W[I-4]);
          K1 = 23.0*F(T, W[I-1][0], W[I-1][1]) - 16.0*F(T-H, W[I-2][0], W[I-2][1]) + 5.0*F(T-2.0*H, W[I-3][0], W[I-3][1]);
          L1 = 23.0*F2(T, W[I-1][0], W[I-1][1]) - 16.0*F2(T-H, W[I-2][0], W[I-2][1]) + 5.0*F2(T-2.0*H, W[I-3][0], W[I-3][1]);
          //W[I] = W[I-1] + H/24.0*K1;
          W[I][0] = W[I-1][0] + H/12.0*K1;
          W[I][1] = W[I-1][1] + H/12.0*K1;

          /* COMPUTE T(I) */
          T = A + I * H;

          /* STEP 7 */
          cout <<"F1 - F2: At time "<< T <<" the solution = "<< W[I][0] << "," << W[I][1] << endl;
          for (int b=0;b<2;b++)
             {
                 str[I-1][b] = W[I][b];
                 myoutput << str[I-1][b]<< " ";
             }
      }

    }
      system("pause");
      return 0;
}
/*  Change function F for a new problem   */
double F(double T, double Y, double Y2)
{
   double f;

   f = -501*Y + 500*Y2 + T*0;
   return f;
}
double F2(double T, double Y, double Y2)
{
   double f;

   f = 500*Y - 501*Y2 + T*0;
   return f;
}