Implement a stiff differential equation solver

[calcmogul]

Sometimes, the physics sim differential equations are stiff enough that the default solver takes way too long because it takes very small steps to stay below the error bound. We should add a higher order solver designed for stiff differential equations.

We should also automatically select the appropriate algorithm for a given problem:
https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Choosing-a-Solver-Algorithm

[calcmogul]

Here's a candidate method: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.Radau.html aka radau5. It's an implicit 5th order method. I haven't found the paper yet, but the scipy docs imply it uses sparse linear algebra, so something bespoke with Eigen::SimplicialLDLT or a problem built in Sleipnir may be warranted.

Here's the Butcher tableau for 5th order Radau IIA: https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Radau_IIA_methods

The gist of implicit methods seems to be rootfinding. For example, solving k₁ = f(yₙ + hk₁) for k₁ (or f(yₙ + hk₁) − k₁) = 0 for k₁) where dy/dt = f(y), then computing yₙ₊₁ = yₙ = hk₁.