I changed the code to the rolling torus example pretty significantly. Little fixes in how the arrows representing the body fixed unit vectors are now fixed. Additionally, I added some code that calculates the kinetic, potential and total energy of the system, and plots it. The scipy solver seems to do a pretty good (not perfect) job of maintaining constant energy. I think tightening up the error tolerances would improve this, but the simulation is well behaved at this point so I didn’t bother.