Use a symplectic integrator.

[idontgetoutmuch]

This preserves the Hamiltonian / Energy but it is O(h) so you may wish to stick with rk4/5 which is O(h^4) for demonstration purposes (I'd recommend pointing out to users that it does not preserve the system). On the other hand, Stan (http://mc-stan.org) uses leapfrog which is O(h^2) and no more difficult to implement than symplectic Euler so maybe it would be worth implementing that for now? Ideally one should write a symplectic order 4 integrator but I don't have time to do that at the moment.

[mstksg]

Thanks for this! I saw the demonstrations in the comment to the blog post and it's pretty awesome. I'm impressed that we can preserve energy pretty much perfectly.

I'm wondering what the best route here is...maybe to offer multiple integrators (as an input enumerable, maybe) so people can decide?

[idontgetoutmuch]

Maybe but it's not a trivial undertaking. I found Ernst Hairer's implementation: https://www.unige.ch/~hairer/software.html and a Julia implementation: https://github.com/DDMGNI/GeometricIntegrators.jl. I certainly would need to refresh my knowledge by ploughing through Iserles' book and also through Hairer's book. It is certainly an interesting project but I need to finish other things before starting on yet something else.