/* compile and link with
gcc -o testcase2 testcase2.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-25, but not for 1e-24 */
  double epsabs = 1e-25;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 1.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 1.0,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}