[Help-gsl] solving ODE in a loop

[Lado Samushia]

I am trying to solve am ode system in a double loop (for different values of
parameters).
This is how my code looks like schematically:
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
int func (double t, const double y[], double f[], void *params){
*********************************************************************
}


int jac (double t, const double y[], double *dfdy, double dfdt[], void
*params){
***********************************************************************
}

int main ( ) {
  double t;
  double h = 1e-6;
  double y[3];
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (3);
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 3);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);

  for ( ) {
    for ( ) {
      t = 1e-6;
      //initialize parameters
      //initialize y[]
      gsl_odeiv_system sys = {func, jac, 3, parameters};
      while ( )
        int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
    }
  }
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);

  return 0;
}
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
It works fine for the first time, but brakes down on second (and consequent)
loop executions (gives "nan" for function values and jumps trough interval
in one step).
Can you tell me how I can fix this problem?
Thank you in advance.