// solve y'(t) = -t*y ,  condition y(0)  = -2
// g++ -std=c++11 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>

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){
   std::cout<<std::fixed<<std::setprecision(16);
   static bool first=true;
   if(first) { 
     std::cout<<std::setw(20)<<"t"<<std::setw(20)<<"gsl solution";
     std::cout<<std::setw(20)<<"exact solution"<<std::setw(20)<<"error\n";
     first = false;
   }
   std::cout<<std::setw(20)<<t<<std::setw(20)<<y<<std::setw(20)<<exact<<std::setw(20)<<y-exact<<"\n";   
}

int main ()
{
  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
  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
  auto 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 (int i = 0; i<n; ++i) {
    auto t = t0+i*dt; // begin of t interval 
    if(gsl_odeiv2_driver_apply (driver, &t, t+dt, y) != GSL_SUCCESS)
      {
	std::cout<<"FAILED near "<<t<<"\n";return 1;
      }      
    exact = -2.0*exp(-0.5*t*t);
    output(t,y[0],exact);
  }  
  gsl_odeiv2_driver_free(driver);  
}