#include <iostream>
#include <boost/numeric/odeint.hpp>
#include <fstream>

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

/* Solves x''(t) = -x(t) - gam*x'(t), split to coupled 
    x'(t) = y(t)
    y'(t) = -x(t) - gam*y(t)
   notation: x(t)=x[0], y(t)=x[1], x'(t) = dxdt[0], y'(t) = dxdt[1] 
*/
void harmonic_oscillator(const state_type &x, state_type &dxdt, const double /*t*/)
{
  const double gam=0.15;
  dxdt[0] = x[1];
  dxdt[1] = -x[0] - gam*x[1];
}


int main()
{
  stepper_type stepper;
  state_type x ={1.0,2.0} ; // initial values, x(0) = 1, x'(0) = 2
  // integrate all the way to final time:
  //integrate_const( stepper , harmonic_oscillator , x , 0.0 , 100.0 , 0.01 );
  //cout<<"final point "<<x<<endl;
  
  // use method do_step to follow the solution step by step
  std::ofstream out("ode.dat");
  const double dt = 0.01;
  for( double t=0.0 ; t<100.0 ; t+= dt ) 
    { 
      stepper.do_step( harmonic_oscillator , x , t , dt );
      std::cout<<t<<" "<<x[0]<<" "<<x[1]<<"\n";
      out<<t<<" "<<x[0]<<" "<<x[1]<<"\n";
    }
  out.close();
}