Design pattern of an ODE class in C++

I've been using different ODE solvers in C++, and I'm familiar with the theory behind standard and more advanced numerical methods for ODEs. What I'd like to understand is what is the "design pattern" of an ODE class. For instance, looking at this question

we notice that the definition of the r.h.s consists of a void function which takes references to z and dzdt as follows:

void ode( const state_type &z , state_type &dzdt , double t ) { 
    dzdt[0] = z[1]; 
    dzdt[1] = -1 * z[0] * w * w; 
} 

and then the integration is carried in the main with

int main() { ... 

    integrate( ode , z , t , 1000 , 0.1 , write_ode ); 
    return 0;

}

Of course such library is really hard-coded, but I just want to grasp the general idea behind the "stepper", let's say for explicit Euler's method.

---

Since the r.h.s is defined like in void ode (...), I imagine that in the stepper part there's a call to void ode(...) that allows to update dzdt. It could be implemented as follows (using the std::vector class)

void do_step(std::vector<double>& z, double tn, double h){
    //tn current time, h time step
    std::vector<double> dzdt(2); 
    ode(tn,y,dzdt);
    z[0] += h*dzdt[0];
    z[1] += h*dzdt[1];
}

i.e.:

• I allocate space for the new solution vector
• compute the new state
• Update z[0],z[1] using Euler's scheme, or whatever.

which can be called with do_step(y,tn,h);

Summarizing, my question is: given the definition of the r.h.s., is this the good way to define a step of a numerical method? Any reference/book with some design techniques of such a problem is highly appreciated