## Examples

The following program solves the second-order nonlinear Van der Pol oscillator equation,

    x"(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0

This can be converted into a first order system suitable for use with the library by introducing a separate variable for the velocity, y = x'(t),

    x' = y
    y' = -x + \mu y (1-x^2)

The program begins by defining functions for these derivatives and their Jacobian,

    #include <stdio.h>
    #include <gsl/gsl_errno.h>
    #include <gsl/gsl_matrix.h>
    #include <gsl/gsl_odeiv.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]*(y[0]*y[0] - 1);
      return GSL_SUCCESS;
    }

    int
    jac (double t, const double y[], double *dfdy,
         double dfdt[], void *params)
    {
      double mu = *(double *)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*mu*y[0]*y[1] - 1.0);
      gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
      dfdt[0] = 0.0;
      dfdt[1] = 0.0;
      return GSL_SUCCESS;
    }

    int
    main (void)
    {
      const gsl_odeiv_step_type * T
        = gsl_odeiv_step_rk8pd;

      gsl_odeiv_step * s
        = gsl_odeiv_step_alloc (T, 2);
      gsl_odeiv_control * c
        = gsl_odeiv_control_y_new (1e-6, 0.0);
      gsl_odeiv_evolve * e
        = gsl_odeiv_evolve_alloc (2);

      double mu = 10;
      gsl_odeiv_system sys = {func, jac, 2, &mu};

      double t = 0.0, t1 = 100.0;
      double h = 1e-6;
      double y[2] = { 1.0, 0.0 };

      while (t < t1)
        {
          int status = gsl_odeiv_evolve_apply (e, c, s,
                                               &sys,
                                               &t, t1,
                                               &h, y);

          if (status != GSL_SUCCESS)
              break;

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

      gsl_odeiv_evolve_free(e);
      gsl_odeiv_control_free(c);
      gsl_odeiv_step_free(s);
      return 0;
    }

The main loop of the program evolves the solution from (y, y') = (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 y.

To obtain the values at regular intervals, rather than the variable spacings chosen by the control function, the main loop can be modified to advance the solution from one point to the next. For example, the following main loop prints the solution at the fixed points t = 0, 1, 2, \dots, 100,

      for (i = 1; i <= 100; i++)
        {
          double ti = i * t1 / 100.0;

          while (t < ti)
            {
              gsl_odeiv_evolve_apply (e, c, s,
                                      &sys,
                                      &t, ti, &h,
                                      y);
            }

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

It is also possible to work with a non-adaptive integrator, using only the stepping function itself. The following program uses the `rk4` fourth-order Runge-Kutta stepping function with a fixed stepsize of 0.01,

    int
    main (void)
    {
      const gsl_odeiv_step_type * T
        = gsl_odeiv_step_rk4;

      gsl_odeiv_step * s
        = gsl_odeiv_step_alloc (T, 2);

      double mu = 10;
      gsl_odeiv_system sys = {func, jac, 2, &mu};

      double t = 0.0, t1 = 100.0;
      double h = 1e-2;
      double y[2] = { 1.0, 0.0 }, y_err[2];
      double dfdy[4], dydt_in[2], dydt_out[2];

      /* initialise dydt_in */
      GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, 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;

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

      gsl_odeiv_step_free(s);
      return 0;
    }

The derivatives and jacobian must be initialised for the starting point t=0 before the first step is taken. Subsequent steps use the output derivatives dydt_out as inputs to the next step by copying their values into dydt_in.