My old, own ODE integration approach fails for three
oscillators model. It probably does so because it has too many
dimensions.
After some deliberation, I decided to kick it out.
Instead, we now use the scipy integration.
However, we need to deal with situations where the system
diverges, where either control or state go to infinity.
We can do this because we use the iterative integration
method. Once an integration step fails, we discard the
previous data and start again, but set the time limit to
something between the largest successful and smallest failing
time.
This way, we should be able to step-by-step approximate the
time limit until which an ill-behaved system remains in an
acceptable state.
The side-effects of the new method are:
1. It may be slower.
2. The number of samples changes dynamically and is no longer
fixed.
3. If we approximate the differential from the data, then the
values we find in the integration result are actually not
real differential values.
They are differences for different interpolation steps of
the function.