/* solves van der pol diffeq

      x'' + mu*x'(x^2 - a^2) - x = 0.

via 4th-order runge-kutta method
also shows how to pass 2 parameters
(in this case, mu & a) to an ode. */


#include <iostream>
#include <cmath>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

using namespace std;


int
     func (double t, const double y[], double f[],
           void *params)
     {

       struct duo{    // this object holds 2 parameters
	 double mu;
	 double a;
       };
       
       duo PARAM;
       PARAM = *(duo *)params;
       //       double mu = *(double *)params;
       f[0] = y[1];                    
       f[1] = -y[0] - PARAM.mu*y[1]*(y[0]*y[0] - pow(PARAM.a, 2.0)); // CHANGE ME
       return GSL_SUCCESS;
     }
     
     int
     jac (double t, const double y[], double *dfdy, 
          double dfdt[], void *params)
     {
       struct duo{
	 double mu;
	 double a;
       };
       
       duo PARAM;
       PARAM = *(duo *)params;


       gsl_matrix_view dfdy_mat 
         = gsl_matrix_view_array (dfdy, 2, 2);
       gsl_matrix * m = &dfdy_mat.matrix; 
       gsl_matrix_set (m, 0, 0, 0.0);
       gsl_matrix_set (m, 0, 1, 1.0);
       gsl_matrix_set (m, 1, 0, -2.0*PARAM.mu*y[0]*y[1] - pow(PARAM.a,2.0));
       gsl_matrix_set (m, 1, 1, - PARAM.mu*(y[0]*y[0] - pow(PARAM.a, 2.0)));
       dfdt[0] = 0.0;
       dfdt[1] = 0.0;
       return GSL_SUCCESS;
     }
     


int main ()
     {
       const gsl_odeiv_step_type * T 
         = gsl_odeiv_step_rk4;
     
       gsl_odeiv_step * s 
         = gsl_odeiv_step_alloc (T, 2);
     
       struct duo{
	 double mu;
	 double a;
       };
       
       duo PARAM;
       PARAM.mu = 0.05;   //CHANGE ME. damping parameter
       PARAM.a = 3.0;     //CHANGE ME. damping parameter


       gsl_odeiv_system sys = {func, jac, 2, &PARAM};
     
       double t = 0.0, t1 = 100.0;   // CHANGE ME. bounds.
       double h = 1e-2;              // CHANGE ME. step size

       // initial conditions: y[0] = position, y[1] = velocity
       double y[2] = { 3.0, 0.0 }, y_err[2];
       double dydt_in[2], dydt_out[2];
     
       /* initialise dydt_in from system parameters */
       GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);
     
       while (t < t1)
         {
           int status = gsl_odeiv_step_apply (s, t, h, 
                                              y, y_err, 
                                              dydt_in, 
                                              dydt_out, 
                                              &sys);
     
           if (status != GSL_SUCCESS)
               break;
     
           dydt_in[0] = dydt_out[0];
           dydt_in[1] = dydt_out[1];
     
           t += h;
	   cout << t << " " << y[0] << " " << y[1] << endl; 

         }
     
       gsl_odeiv_step_free (s);
       return 0;
     }