/*
 * The program solves the second-order nonlinear 
 * Van der Pol oscillator equation,
 *
 *   u''(t) + \mu u'(t) (u(t)^2 - 1) + u(t) = 0
 *
 * This can be converted into a first order system 
 * by introducing a separate variable for the velocity, 
 * v = u'(t),
 *
 *   u' = v
 *   v' = -u + \mu v (1-u^2)
 *
 * The program begins by defining functions for these 
 * derivatives. The main function uses 
 * driver level functions to solve the problem. The program 
 * evolves the solution from (u, v) = (1, 0) at t=0 to t=100. 
 * The step-size h is automatically adjusted by the controller 
 * to maintain an absolute accuracy of 10^{-6} in the function 
 * values (u, v). 
 */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int
func (double t, const double y[], double f[], void *params)
{
  double mu = *(double *) params;
  f[0] = y[1];
  f[1] = -y[0] + mu*y[1]*(1. - y[0]*y[0]);
  return GSL_SUCCESS;
}

int
main (void)
{
  size_t neqs = 2;          /* number of equations */
  double eps_abs = 1.e-6, eps_rel = 0.;  /* desired precision */
  double stepsize = 1e-6;   /* initial integration step size */
  double mu = 10.;          /* the parameter of the oscillator */
  double t = 0., t1 = 100.; /* time interval of the evolution */
  /* 
   * Initial conditions 
   */
  double y[2] = { 1., 0. };

  /*
   * Allocate integrator
   *
   * We use explicit embedded Runge-Kutta-Fehlberg (4,5) method. 
   * This method is a good general-purpose integrator.
   */
  gsl_odeiv2_step    *s = gsl_odeiv2_step_alloc (gsl_odeiv2_step_rkf45, neqs);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
  gsl_odeiv2_evolve  *e = gsl_odeiv2_evolve_alloc (neqs);
    
  gsl_odeiv2_system sys = {func, NULL, neqs, &mu};
    
  /*
   * Evolution loop 
   */
  while (t < t1) {
      int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &stepsize, y);

      if (status != GSL_SUCCESS)
          break;

      printf ("% .5e % .5e % .5e\n", t, y[0], y[1]);
  }

  gsl_odeiv2_evolve_free (e);
  gsl_odeiv2_control_free (c);
  gsl_odeiv2_step_free (s);

  return 0;
}