Currently, opinf.models.ContinuousModel.predict() wraps scipy.integrate.solve_ivp(). It would be nice to be able to pass a custom time-stepper, probably as the method attribute. This is also important for certain types of models, like Hamiltonian systems, which require symplectic integrators.

• Write an IntegratorTemplate somewhere (new opinf.integrate submodule?)
• Update predict() so that method can be an integrator object.
• Write a few common integrators as examples (forward/backward Euler, IMEX, etc.).

Need to think about how the integrator should interact with the list of operators or the model's rhs() method.