// solve y'(t) = -t*y ,  condition y(0)  = -2
// compile :
// g++ --std=c++0x gsl_ode_simple.cpp `gsl-config --libs`
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
using namespace std;
int func (double t, const double y[], double f[], void *params){
  f[0] = -t*y[0] ; // y[0] = y,  f[0] = y' = dy/dt = -t*y 
  return GSL_SUCCESS;
}
void output(double t, double y, double exact){
   cout<<fixed<<setprecision(16);
   static bool first=true;
   if(first) { 
     cout<<setw(20)<<"t"<<setw(20)<<"gsl solution";
     cout<<setw(20)<<"exact solution"<<setw(20)<<"error\n";
     first = false;
   }
   cout<<setw(20)<<t<<setw(20)<<y<<setw(20)<<exact<<setw(20)<<y-exact<<endl;   
}

int main (void)
{
  double exact;
  const int n=50;                 // # of points
  double t0 = 0.0, t1 = 10.0;     // time start and end
  double dt=(t1-t0)/(n-1);        // time step
  vector<double> tim(n-1);        // time interval start points
  {
  int i=0;
  for( auto & t: tim) {t=t0+i*dt; ++i;}  // t0,t0+dt,..,t1-dt
  }
  double y[1] = {-2.0};           // initial value; table with one entry

  gsl_odeiv2_system sys = {func, NULL, 1, NULL}; // not using a Jacobian, hence NULL pointer
  gsl_odeiv2_driver* driver =  
    gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-13, 1e-13, 0.0);
  
  output(t0,y[0],y[0]);
  for (auto t: tim) {
    int status = gsl_odeiv2_driver_apply (driver, &t, t+dt, y);      
    if (status != GSL_SUCCESS) {cout<<"FAILED near "<<t<<endl;return 1;}
    exact = -2.0*exp(-0.5*t*t);
    output(t,y[0],exact);
  }  
  gsl_odeiv2_driver_free(driver);  
  return 0;
}