Generic Dosing Mechanism

For many practical pharmacometric models the integrate_ode function is pretty limiting as it only solves the ODE system at given times but cannot modify the system at those times. A more generic function would require events, which pair an event time and some used-defined function that modifies the state of the ODE at that time, for example

void event_action(real[] states, real[] params, real[] x).

These events would then be collected into a dynamically-sized vector of events,

event_table events;
event_table.add(1, increment_1);
...
event_table.add(493, multiply_2);

which would be passed into a new integrate function,

integrate_ode_with_events(ode, y0, t0, ts, events, theta, x_r, x_i).

This new function then concatenates the vector of observation times and the
event times, removes duplicates, and orders the unique times. The solver is
run to each time, with the resulting state saved if an observation if necessary
and then input to the event functor if necessary.

Thoughts? This would require introducing a type with a method which would
be new to the language. Maybe just having an add_event function with an
event_table and function as arguments?

Start by locating the existing integrate_ode implementation and the surrounding ODE interface. Compare its handling of observation times with the proposed event_table and integrate_ode_with_events API; done would require a settled event representation and behavior for ordered, deduplicated observation and event times.