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

Mat_DP *yy_p;

void derivhom(const DP x, Vec_I_DP &y, Vec_O_DP &dydx);
void derivinh(const DP x, Vec_I_DP &y, Vec_O_DP &dydx);
double yexact(const double x);

#define NN 50
int nvar;
const DP ee = exp(1.0);
DP x, x1, x2, h, v1, v2, kappa;

int main(int nbargs, char* args[])
{
  if (nbargs==1) {
     // No argument provided, here ignored 
     nvar = 2;
     x1 = 0.0;   // Starting point
     x2 = 1.0;   // End Point

     v1 =  0.0;
     v2 =  1.0;

     Vec_IO_DP yy1(nvar), yy2(nvar), yy(nvar), f(nvar);
     Vec_DP v(nvar);

     v[0] = 1.0;  // b.v. LHS 
     v[1] = v1;   // Initial guess

     bool check;
     // NR::newt(v,check,NR::shoot);
     // if(check) cout << "Convergence problem " << endl;
     // cout << "v[0] = " << v[0] << ",  v[1] = " << v[1] << endl;
     // NR::odeint(yy,x1,x2,EPS,h,h,nok,nbad,derivinh,NR::rkqs);
     x = x1; 
     h = (x2 - x1)/NN;
     yy1[0] = v[0];  // i.c. for the original problem
     yy1[1] = v1; 
     yy2[0] = 0.0;   // i.c. for the homogeneous problem
     yy2[1] = v2;  

     yy_p = new Mat_DP(3*nvar,NN+1);
     Mat_DP &yp = *yy_p; 

     yp[0][0] = yy1[0];
     yp[1][0] = yy2[0];

     for(int j=1;j<=NN;j++) {
        x += h;
        derivinh(x,yy1,f);
        NR::rk4(yy1, f, x, h, yy1, derivinh);
        derivhom(x,yy2,f);
        NR::rk4(yy2, f, x, h, yy2, derivhom);
        yp[0][j] = yy1[0];
        yp[1][j] = yy2[0];
     }
     // linear combination: compute the constant from the last values
     kappa = (1.0 - yp[0][NN])/yp[1][NN]; 
     x = x1;
     cout << "#  x  " << setw(10) << " y_appr" << setw(16) << " y_ex";
                           cout << setw(20) << " abs(y_a - y_e)" << endl;  
     for(int j=0;j<=NN;j++) {
        yp[2][j] = yp[0][j] + kappa*yp[1][j];
        cout << x << setw(16) << yp[2][j] << setw(16) << yexact(x);
           cout << setw(16) << abs(yp[2][j] - yexact(x)) << endl;
        x += h;
      }
  }
  return 0;
}      

void derivhom(const DP x, Vec_I_DP &y, Vec_O_DP &dydx)
{
  dydx[0] = y[1];
  dydx[1] = exp(-x)*y[1] + y[0];
}

void derivinh(const DP x, Vec_I_DP &y, Vec_O_DP &dydx)
{
  dydx[0] = y[1];
  dydx[1] = exp(-x)*y[1] + y[0] + 1.0;
}

double yexact(const double x)
{ 
  const double c1 = (1.0 + ee - 2.0*ee*ee/exp(1.0/ee) )/(1-ee);
  const double c2 = 2.0*ee;

  return( c1*(1.0-exp(x)) + c2*exp(x-exp(-x)) - exp(x) );
}

/*
 * GNUPLOT Output:
 * plot "grxx113.dat" u 1:2 w l, "grxx113.dat" u 1:3 w l # shoot sol., ex. sol
 * plot "grxx113.dat" u 1:4 w l                          # abs. error
 */