C++ class member function to GSL ODE solver

I am struggling with -- apparently -- a well known problem. I have an ODE system which is defined by a class member function, and I want to solve/integrate it by one of the GSL solvers. I.e., say my model is defined by:

class my_model{
...
public:
    int NEQ = 4;
    double y[4], dydt[4];
    double params[25]; 
    int ode_gsl(double t, const double y[], double dydt[], void * params);
    ...
};

int my_model::int ode_gsl(double t, const double y[], double dydt[], void * params){
    dydt[0] = params[1]*params[0]*y[0] - y[1];
    dydt[1] = ...;
    ...
    return GSL_SUCCESS;    
}

Then in my integration routine:

int main(){
    my_model chemosc;

   // Parameters for the Integrator
   double HSTART = 1e-3;
   double ATOL = 1e-8;
   double RTOL = 1e-6;

   // Instantiate GSL solver
   gsl_odeiv2_system sys = {&chemosc.ode_gsl, nullptr, chemosc.NEQ, chemosc.params};
   gsl_odeiv2_driver * d;
   d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, HSTART, ATOL, RTOL);

   double twin[2] = {0.,60.};
   double dt = 1e-3;
   double t = twin[0], t1 = twin[1];
   long int NSTEP = (long int)((t1-t)/dt)+1; // +1 if you start counting from zero...
   int NEQ = 4;
   long int NUMEL = (NEQ+1)*NSTEP; // number of elements for solution

   int i = 0,j;
   do{
      double ti = t + dt; // t is automatically updated by the driver
      printf("\n%.3f\t%.3f\t%.3f t%.3f",astro.y[0],astro.y[1],astro.y[2],astro.y[3]);
      int status = gsl_odeiv2_driver_apply (d, &t, ti, astro.y);
      ...
      }
...
}

Compiling the above code creates the somehow well known error that GSL requires a pointer-to-function whereas I am passing a pointer-to-member function, i.e.:

error: cannot convert ‘int (chemosc::*)(double, const double*, double*, void*)’ to ‘int (*)(double, const double*, double*, void*)’

I found several questions related to this topic: Q1, Q2, Q3, Q4, Q5, Q6, but seriously, the answer there are hard to follow. Declaring as static my member function, has the downside that the compiler requires me to also declare as static all the member parameters. Using static_cast as suggested here, result in all bunch of segmentation errors (but I assume that I did something wrong in the implementation, inasmuch as the directions in that post are quite cryptic). It would be nice to have a once-for-all clear and as simple as possible working solution for this problem (possibly without using boost libraries).

Something like this:

class my_model
{
    private: int
    ode_gsl_impl(double const t, double const * const y, double * const dydt);

    public: static int
    ode_gsl(double const t, double const * const y, double * const dydt, void * const opaque)
    {
        assert(opaque);
        return(static_cast<my_model *>(opaque)->ode_gsl_impl(t, y, dydt));
    }
};

gsl_odeiv2_system sys
{
    &my_model::ode_gsl
,   nullptr
,   chemosc.NEQ
,   reinterpret_cast<void *>(::std::addressof(chemosc))
};

I would like to mention that your original code suffers from name collisions between callback parameter names and class field names.