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

// The function calculating the rhs of the ODE

int
func (double t, const double y[], double dydt[],
      void *params)
{

  // params is casted to 'double *' and assigned to my_array.
  // Now 'my_array' is identical to the variable 'array' in main ()
  // as intended.
  
  const double * my_array = (const double *) params;

  // Evaluate the function
  
  dydt [0] = t + my_array [0] + my_array [1] + my_array [2];
  return GSL_SUCCESS;
}
     
int
main (void)
{
  // This is the array that is used as a parameter of the function
  // that calculates the rhs of the ODE
  
  double array [3] = {2, 4, 6};

  // Initialise ODE solver

  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rkf45;
     
  gsl_odeiv_step * s 
    = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c 
    = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e 
    = gsl_odeiv_evolve_alloc (1);
     
  gsl_odeiv_system sys = {func, 0, 1, array};
     
  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y [1] = { 0.0 };
     
  while (t < t1)
    {
      int status = gsl_odeiv_evolve_apply (e, c, s,
					   &sys, 
					   &t, t1,
					   &h, y);
     
      if (status != GSL_SUCCESS)
	break;
    }
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}