gsl_function alternative for c++

I am switching from C to C++, and I would like to optimally use the additional features that become available and avoid things considered 'C-style' like void * pointers. Specifically, I am trying to make a gsl_function-like interface (NOT a wrapper to use gsl in C++).

In C, I wrote several routines for root-finding, integration ... that use a gsl_function-like interface to pass mathematical functions to these routines. This interface looks like this:

struct Function_struct
{
  double (* p_func) (double x, void * p_params);
  void * p_params;
};
typedef struct Function_struct Function;

#define FN_EVAL(p_F,x) (*((p_F)->p_func))(x,(p_F)->p_params)

and can be used in the following way:

struct FuncParams_struct { double a; double b; double c; };

double my_sqrt(double x, void * p) {
    struct FuncParams_struct * p_params = (struct FuncParams_struct *) p;
    double a = p_params->a;
    double b = p_params->b;
    double c = p_params->c;

    return a/sqrt(x+b)+c;
}

Function My_Sqrt;
My_Sqrt.p_func = &my_sqrt;
struct FuncParams_struct my_params = { 1.0, 1.0, 0.0 };
My_Sqrt.p_params = &my_params;

// Call the function at a certain x
double result_at_3 = FN_EVAL(&My_Sqrt, 3);

// Pass My_Sqrt to an integration routine 
// (which does not care about the parameters a,b,c 
// and which can take any other Function to integrate)
// It uses FN_EVAL to evaluate MySqrt at a certain x.
double integral = integrate(&My_Sqrt);

// Easily change some of the parameters
My_Sqrt.p_params->a=3.0;

As you can see, this allows me to create an instance of the Function structure, which contains a pointer to some function parameters (which are typically contained within another structure) and a pointer to a 'normal' C function. The great thing about this is that the routines that use Function only need to know that it is a function that takes a double and returns a double (they use the FN_EVAL macro). They do not have to care about the type and number of parameters. On the other hand, I can easily change the values stored in the parameters from outside of the routines.

I would now like to have the features highlighted above in a C++ Function. I searched a lot for what would be the best way to get this, but I could not find it (partly because I am a bit overwhelmed by the possibilities of C++ as compared to C). Thus far, I was thinking to make a class Function. This class could than store the actual function definition and it could be made callable by defining the operator(), such that again routines do not have to care about the parameters. The operator() could thus be used instead of the FN_EVAL macro.

The things that I could not yet decide/find are:

• how I would store the parameters. I was thinking about using a template typename. However, as far as I understand, then also the class itself needs to be a template, in which case also the routines should accept a templated class. I do not want to use a void *.
• how I would change the parameters that are stored in my class Function. Should I make them public such that I can easily change them, or should I keep them private and write some interface to access them? In the latter case, how should I go about this? Since I will be considering a large variety of parameters (both in number and in type).

Things to consider:

• I found a lot of questions and answers about wrappers to use gsl-routines within C++. That is not what I am looking for.
• I also looked at the STL <functional>. As far as I could figure out, it does not fulfill my requirements for the parameters that I can store within the Function. However, I can be wrong about this.
• My functions can be rather complex (i.e., not just one-line things like the example above) and the parameters can themselves contain a mixture of doubles, ints, structures...

• I would like to have the possibility to extend the Function class to higher dimensions, i.e. f(x,y,z). In C, I had for example

  struct FunctionThree_struct
  {
    double (* p_func) (double x, double y, double z, void * p_params);
    void * p_params;
  };
  typedef struct FunctionThree_struct FunctionThree;
  #define FN_THREE_EVAL(p_F,x,y,z) (*((p_F)->p_func))(x,y,z,(p_F)->p_params)
  
  

• I do care about efficient function calls. Especially for higher dimensional integrals, the functions will be called millions of times.