Some of that Python is strange for numerical Python practices, yes, but I was aiming for close floating-point output equivalence. Rust ndarray provides the standard set of matrix ops, but I keep getting different results (\propto 1e-5 on a stable basis, which isn't quite a proof of incorrectness, but....) when I use too much numpy stuff like linalg.norm, etc. (ndarray slices in Rust are annoying to use because of lifetimes and such.)

I did my dissertation on symplectic exponential Runge-Kutta schemes with stable numerics for the big matrix exponentials needed (which have this specific structure that allow for theorems to be proven). But I didn't have time to write any code at all. I wonder if open source ODE solvers are getting good high-order symplectic methods by now...