#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
using namespace std;


class F{
  public:
    F(){;};
    static int f(double t, const double* y, double* yprime, void* s ){
      return 0;
    }
};


int main(){

  double t = 0;
  double nextt = 100;
  double internalStepSize_ = 0.01;
  const gsl_odeiv_step_type* gslStepType_ = gsl_odeiv_step_rkf45;
  gsl_odeiv_step* gslStep_ = 0; 
  gsl_odeiv_control* gslControl_;
  gsl_odeiv_evolve* gslEvolve_;
  gsl_odeiv_system gslSys_;
  double *y_ = 0;
  double absAccuracy_ = 1.0e-6;
  double relAccuracy_ = 1.0e-6;
  int nVarMols_ = 19;
  F *g = new F();
  
  y_ = new double[ nVarMols_ ];
  double d;
  memcpy( y_, &d, nVarMols_ * sizeof( double ) );
  
  if ( gslStep_ ) gsl_odeiv_step_free( gslStep_ );
  gslStep_ = gsl_odeiv_step_alloc( gslStepType_, nVarMols_ );
  gslControl_ = gsl_odeiv_control_y_new( absAccuracy_, relAccuracy_ );
  gslEvolve_ = gsl_odeiv_evolve_alloc( nVarMols_ );
  gslSys_.function = &F::f;
  gslSys_.jacobian = 0;
  gslSys_.dimension = nVarMols_;
  gslSys_.params = g;
  
  if ( gslStep_ ) gsl_odeiv_step_reset( gslStep_ );
  
  int status = gsl_odeiv_evolve_apply ( gslEvolve_, gslControl_, gslStep_, &gslSys_, &t, nextt,	&internalStepSize_, y_);
  cout << *y_ << endl;
  return 0;
}