#include <cstdlib>
#include <ctime>
   // for time measuring
#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;

int jj = 0;              // (global) counter for RHS computations
const double a = 0.2;    // global parameter
const double b = 0.03;
const double c = 0.1;
const double d = 0.03;
const double y10 = 6.0;  // initial condition
const double y11 = 15.0;
const double TE = 70.0;  // end point
const DP TOL = 1.0e-6; 

void deriv(const DP t, Vec_I_DP &y, Vec_O_DP &f)
  // RHS of ODE
{
  f[0] = -a*y[0] + b*y[0]*y[1];
  f[1] =  c*y[1] - d*y[0]*y[1];
  jj++;             // counter for RHS calls
}

void simpson( int k, Vec_I_DP &y, const DP t, const DP h, const DP TOL, 
                     Mat_IO_DP yp, Vec_IO_DP xp, Vec_O_DP &yout )
{
  int i, j, nmax = 4;
  double delta;

  int nvar=y.size();
  Vec_DP u(nvar), v(nvar), w(nvar);
  Vec_DP ff(nvar);

  // constant part for the current step 
  for (i=0; i<nvar; i++) {
     w[i] = yp[i][k-2] + (yp[nvar+i][k-2] + 4.0*yp[nvar+i][k-1])*h/3.0; 
     v[i] = yp[i][k-1];
  }
  delta = 1.0e3;
  j = 0;             // iteration counter 
     
  while ( delta > TOL && j < nmax )  // iteration loop
  {
     deriv(t, v, ff);
     delta = 0.0;
     for(i=0;i<nvar;i++) {
        u[i] = w[i] + ff[i]*h/3.0;
        delta += fabs(u[i] - v[i]);
     }
     for(i=0;i<nvar;i++) v[i] = u[i];
     j++;
  }
  for(i=0;i<nvar;i++) {
     yout[i] = u[i];
     yp[i][k] = u[i];
  }
  deriv(t, u, ff);
  for(i=0;i<nvar;i++) yp[nvar+i][k] = ff[i];
  // cout<<" k = "<< k <<",  #it = "<< j <<",   delta = "<< delta << endl;
  return;
}

// ======================================================================

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

   hh = TE/(double)NN;
   tt = 0.0;
   y[0] = y10; 
   y[1] = y11;

   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: " << endl;
   //    cout << "# j = 0 " << "\t" << xp[0] << "\t" << yp[0][0] << "\t";
   //      cout << yp[0][1] << "\t" << yp[N][0] << "\t" << yp[N][1] << endl; 
 
   for (j=1; j<2; j++) {   // startig procedure (Runge-Kutta)
      deriv(tt,y,f);
      NR::rk4(y, f, tt, hh, yout, deriv);
      tt += hh;
      xp[j] = tt;
      deriv(tt,yout,f); // function value of solution point
      for (i=0;i<N;i++) {
          yp[i][j]=yout[i];
          yp[N+i][j]=f[i];
      }
   }
   // cout << "# -> end of starting procedure ..." << endl;

   // multistep method ( Simpson rule without predictor )
   double cpuTotal = 0.0;
   clock_t startTime = clock();

   for (j=2; 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: Simpson 
      simpson(j, y, tt, hh, TOL, yp, xp, yout);
                        
      tt += hh;
      xp[j] = tt;
      deriv(tt,yout,f);
      for (i=0;i<N;i++) {
          y[i] = yout[i];
          yp[i][j]=yout[i];
          yp[N+i][j]=f[i];
      }
   }
   clock_t endTime = clock();
   double cpu = double(endTime - startTime) / CLOCKS_PER_SEC;

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

   cout << endl << "#-> numerical solution:" << endl;
   cout << "#  j  |   t     |     y1(t)           y2(t)    | ";
      cout << " f(t,y1(t)     f(t,y2(t)  " << endl;
   cout << "# "; for (i=0;i<72;i++) cout << "="; 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] << setw(14) << yp[2][j] << setw(14) << yp[3][j] << endl;
   }

   cout << endl << " # number of RHS computations: " << jj << endl;
   cout << " # cpu time: " << cpu << " [sec.]" << endl;
   // cout << "CLOCKS_PER_SEC " << CLOCKS_PER_SEC << endl;

   delete xp_p, yp_p;
}