The present time-stepping scheme is the 4th order Gauss-Legendre. The 1e-10 tolerance is not the desired accuracy of the solution but used for the Newton iterations which solve the implicit scheme to machine precision in order to exactly preserve the energy of the conservative dynamical system being integrated.

The goal is to illustrate how one can obtain an engineering tolerance approximation that still preserves energy without needing to compute a high precision solution. In particular, all approximations in the time resolution study preserve energy to within rounding error even the inaccurate ones.

However, I’m admittedly a beginner with IRK methods.

I have looked through your code for the 8-stage IRK scheme. Since I’m integrating a system of three ODE’s my understanding is that this would lead to 24 variables that need to be implicitly solved to machine precision at each time step. The system isn’t particularly stiff, but the idea is to relax accuracy requirements as much as possible while continuing to exactly preserve energy.

I’ll see what I can do with the code you linked. These mathematical pointers and ideas are greatly appreciated. Thanks!