// solve y'(t) = -t*y ,  condition y(0)  = -2
// Adaptive integration using Boost::numeric::odeint
// compile :
// g++ --std=c++11 boost_ode_simple.cpp
#include <iostream>
#include <iomanip>
#include <vector>
#include <boost/numeric/odeint.hpp>

using state_type= std::vector<double>;
using stepper_type= boost::numeric::odeint::runge_kutta_cash_karp54<state_type>;

namespace my{
  using std::cout;
  void system(const state_type& y, state_type& dydt, const double t)
  {
    dydt[0] = -t*y[0];
  }
  void output(double t, double y, double exact){
    cout<<std::fixed<<std::setprecision(16);
    static bool first=true;
    if(first) { 
      cout<<std::setw(20)<<"t"<<std::setw(20)<<"Boost solution";
      cout<<std::setw(20)<<"exact solution"<<std::setw(20)<<"error\n";
      first = false;
    }
    cout<<std::setw(20)<<t<<std::setw(20)<<y<<std::setw(20)<<exact<<std::setw(20)<<y-exact<<"\n";   
  }
}

int main()
{  
  const double e_abs=1e-15,e_rel=1e-15; // absolute and relative error goal
  auto stepper = boost::numeric::odeint::make_controlled<stepper_type>(e_abs, e_rel) ;
  state_type y ={-2.0} ; //condition y(0)=-2 gives y(t) = -2.0*exp(-0.5*t*t)
  const int n=50; // solve 50 points
  const double t1 = 0.0, t2=10.0;
  const double dt = (t2-t1)/(n-1);

  for(int i = 0; i<n ; ++i) 
    { 
      auto t = t1+i*dt;      
      auto exact = -2.0*exp(-0.5*t*t);
      my::output(t,y[0],exact);
      auto steps =integrate_adaptive(stepper, my::system, y, t, t+dt, 0.01);
      // cout<<"steps "<<steps<<endl; // how many sub-division steps were needed
    }
}