Heun 2nd order coupling

[maedoc]

Describe the bug

In the current integrators, the coupling has been evaluated only once per step, while the Heun method should evaluate it twice. For many models, this isn't an issue, but some models have highly nonlinear coupling functions (e.g. Kuramoto) and convergence suffers significantly compared to custom code (10x smaller time step required).

Steps to reproduce

Any model even with linear coupling, e.g. Kuramoto or MPR.

Expected results

Faster convergence.

Actual results

Slow convergence.

Additional information

This will require a change to data flow as the integrator does not currently receive the coupling function instance as an argument.