#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cmath>
#include "nr.h"
using namespace std;

#define PI 3.141592653589793

Vec_DP *xp_p;
Mat_DP *yp_p;

const int n_ord = 3;
const double a = 8.0;   // global parameter
const double ya = 2.0;  // initial condition

double yex(double t)
{
  // analytical solution of the ODE 
  double a2;
  a2 = 1.0/(1.0+a*a);
  return( (ya - a2)*exp(-t) + a2*( cos(a*t) + a*sin(a*t) ) );
}

void deriv(const DP t, Vec_I_DP &y, Vec_O_DP &f)
  // RHS of ODE
{
  f[0] =  cos(a*t) - y[0];
}

void adamsbash(int k, Vec_IO_DP &y, const DP t, const DP h, 
              Mat_IO_DP yp, Vec_IO_DP xp, const int n_ord, Vec_O_DP &yout )
{
  const int NMAXORD = 5;
  static DP beta[NMAXORD][NMAXORD] = {{1.0, 0.0, 0.0, 0.0, 0.0},
              {(3.0/2.0), (-1.0/2.0), 0.0, 0.0, 0.0},
              {(23.0/12.0),(-16.0/12.0),(5.0/12.0), 0.0, 0.0},
              {(55.0/24.0),(-59.0/24.0),(37.0/24.0),(-9.0/24.0), 0.0},
              {(1901.0/720.0),(-2774.0/720.0),(2616.0/720.0),(-1274.0/720.0),
                                                            (251.0/720.0)}};
  int i,j;

  int nvar=y.size();
  Vec_DP sum(nvar);

  if (n_ord > NMAXORD || n_ord < 1) 
                  NR::nrerror("Formula not present");

  for (i=0; i<nvar; i++) sum[i] = 0.0;

  for (j=0; j < n_ord; j++) {
      for (i=0; i<nvar; i++) 
           sum[i] += beta[n_ord-1][j]*yp[nvar+i][k-j-1];
  }
  for (i=0; i<nvar; i++) {
          yout[i] = y[i] + h*sum[i];
          yp[i][j+1] = yout[i];
  }
}

int main()
{
   const int N=1;       // dimension of system
   int i, j, NN =100;   // ## plot points
   DP hh, tt;
   Vec_DP y(N), yout(N), yy(N), f(N);

   hh = 2.0*PI/(double)NN;
   tt = 0.0;
   y[0] = ya; 

   xp_p = new Vec_DP(NN);
   yp_p = new Mat_DP(2*N,NN);
   Vec_DP &xp = *xp_p;
   Mat_DP &yp = *yp_p;

   // starting procedure
   xp[0] = tt; 
   deriv(tt,y,f);
   for (i=0;i<N;i++) {
      yp[i][0] = y[i];
      yp[N+i][0] = f[i];
   }
   cout << "Start: " << "n_ord = " << n_ord << endl;
       cout << "j  = 0 " << "\t" << xp[0] << "\t" << yp[0][0];
                                    cout << "\t" << yp[N+0][0] << endl; 
 
   for (j=1; j < n_ord; j++) {
      deriv(tt,y,f);
      NR::rk4(y, f, tt, hh, y, deriv);
      tt += hh;
      xp[j] = tt;
      deriv(tt,y,f); // function value of solution point
      for (i=0;i<N;i++) {
          yp[i][j]=y[i];
          yp[N+i][j]=f[i];
      }
  
   }
   // cout << " -> end of starting procedure ..." << endl;

   // multistep method ( Adams-Bashforth, order <n_ord> )
   for (j=n_ord; j<NN; j++) {
      tt = xp[j-1];
      for (i=0;i<N;i++) {
         y[i] = yp[i][j-1];
         f[i] = yp[N+i][j-1];
      }
      // multistep method: Adams-Bashforth
      adamsbash(j, y, tt, hh, yp, xp, n_ord, yout);
                        
      tt += hh;
      xp[j] = tt;
      deriv(tt,yout,f);
      for (i=0;i<N;i++) {
          yp[i][j]=yout[i];
          yp[N+i][j]=f[i];
      }
      // cout << "j = " << j << "\t" << xp[j] << "\t"
      //                           << yp[0][j] << "\t" << yp[N+0][j] << endl;
   }

// end of the main procedure
// Simple output for GNUPLOT:

   cout << "#-> analytical solution: " << endl;
   cout << "#  j |   t     |   yex(t)      |  f(t,yex(t)   " << endl;
   cout << "#======================================================= " << endl;
   tt = 0.0;    
   for (int j=0;j<NN;j++) {
     yy[0] = yex(tt);
     deriv(tt, yy, f);
     cout << j << setw(12) <<  tt << setw(14) << yy[0] << setw(14);
                                                  cout << f[0] << endl;
     tt += hh;
   }

   cout << endl << "#-> numerical solution;   n_ord = " << n_ord << endl;
   cout << "#  j |   t     |   yp(t)      |  f(t,yp(t)       " << endl;
   cout << "# ======================================================= " << endl;
   for (int j=0;j<NN;j++) {
     cout << j << setw(12) << xp[j] << setw(14) << yp[0][j] << setw(14);
                                                  cout << yp[1][j] << endl;
   }

   cout << endl << "#-> error:   n_ord = " << n_ord << endl;
   cout << "#  j |   t    |   yp(t)     |  yex(t)    | yp(t) -  yex(t)";
   cout << endl;
   cout << "# ============= =========================================" << endl;
   for (int j=0;j<NN;j++) {
      yy[0] = yex(xp[j]);
      cout << j << setw(12) << xp[j] << setw(14) << yp[0][j] << setw(14);
                 cout << yy[0] << setw(15) << abs(yp[0][j]-yy[0]) << endl;
   }
   delete xp_p, yp_p;
}