/*
 * The program solves the first-order nonlinear equation,
 *
 *   y'(t) = u(t)^2  - u(t)^3
 *
 *   y[0] = 0.0001, 0 <= t <= 20000
 */

#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)
{
  f[0] = y[0]*y[0]*(1. - y[0]);
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 1, 1);
  gsl_matrix *m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, y[0]*(2. - 3*y[0]) );
  dfdt[0] = 0.0;
  return GSL_SUCCESS;
}

int
main (void)
{
  /*
   * Explicit embedded Runge-Kutta-Fehlberg (4,5) method. 
   * This method is a good general-purpose integrator.
   */
  /*const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rkf45;*/
  
  /*
   * Implicit Bulirsch-Stoer method of Bader and Deuflhard. 
   * The method is generally suitable for stiff problems. 
   * This stepper requires the Jacobian.
   */
  const gsl_odeiv2_step_type *T = gsl_odeiv2_step_bsimp;
    
  gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(T, 1);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(1e-9, 0.0);
  gsl_odeiv2_evolve *e = gsl_odeiv2_evolve_alloc(1);

  gsl_odeiv2_system sys = {func, jac, 1, NULL};

  double t = 0., t1 = 20000.;
  double h = 1e-6;
  double y[1] = { 0.0001 };

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

      if (status != GSL_SUCCESS) {
          printf ("Troubles at % .5e % .5e \n", t, y[0]);
          break;
      }

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

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

  return 0;
}