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

// Number of steps for the Euler, mod. Euler, Heun2,and RK4 method
#define M 100

const double x0 = 3.5;    // initial condition
const double a = 4.1;     // parameter
const double pi = 4.0*atan(1.0);
const double T = 2.0*pi;  // Length of the intervall

// diff. equation
void yprime(const DP t, Vec_I_DP &y, Vec_O_DP &f)
{
  f[0] = y[0]*a*sin(a*t);
}

// analytic solution
double solution(double tau)
{
  double sol;

  sol =  x0*exp(1.0-cos(a*tau));
  return sol;
}


int main(int argc, char *argv[])
{
  const int N=1;
  double yex;
  DP  h, t;
  Vec_DP y1(N), y2(N), y3(N), y4(N), f1(N), f2(N), f3(N), f4(N), yerr(N);

  t = 0.0;
  // step size
  h = T/M;
  // initial condition for different methods
  y1[0] = x0;
  y2[0] = x0;
  y3[0] = x0;
  y4[0] = x0;

  cout << "#  " << "   t   " << "     ex_sol.  " << "   Euler   ";
       cout << "   mod_Euler  " << "  Heun_(2)  " << "  RK_(4)  " << endl;
  std::cout.precision(4);
       cout << setw(6) << t << setw(16);
       std::cout.precision(8);
  cout << x0 << setw(12) << y1[0] << setw(12) << y2[0] << setw(12);
                            cout << y3[0] << setw(12) << y4[0] << endl;
 
  for(int j=0;j<M;j++) {
     yprime(t,y1,f1);
     NR::euler(y1,f1,t,h,y1,yprime);

     yprime(t,y2,f2);
     NR::eulermod(y2,f2,t,h,y2,yprime);

     yprime(t,y3,f3);
     NR::heun2(y3,f3,t,h,y3,yprime);

     yprime(t,y4,f4);
     NR::rk4(y4,f4,t,h,y4,yprime);
 
     t+=h;
     yex = solution(t);
     std::cout.precision(4);
     cout << setw(6) << t << setw(16);
     std::cout.precision(8);
     cout << yex << setw(12) << y1[0] << setw(12) << y2[0] << setw(12);
                            cout << y3[0] << setw(12) << y4[0] << endl; 
  }

}
// Output of the result: 'grxx044 > gr044.dat'
//
// Graphics: gnuplot 
// Solution:
//     plot "gr044.dat" u 1:2 w l	
// Error ( e.g. for modEuler and for Heun(2) )
//     plot "gr044.dat" u 1:(abs($4-$2)), "gr044.dat" u 1:(abs(($5-$2))
//