dirk methods for higher order time integration

Instead of relying on time centering for higher order time integration, which can be numerically unstable and induce wringing at times, perhaps we can implement an optional higher order time integration method.

This would involve evaluating the jacobian matrix atleast one additional time per solver iteration, but could potentially enhance solver convergence in certain scenarios. It remains to be seen.

Ideally, we would adopt some sort of diagonally implicit Runga Kutta (DIRK) method? Perhaps a simple two stage DIRK method is enough, and can offer unconditional stability similar to the default backward Euler MESA currently uses.

If we decide to tackle this, these might be some helpful resources:

General review: section 4.1 in https://ntrs.nasa.gov/api/citations/20160005923/downloads/20160005923.pdf for (ESDIRK2 and SDIRK2)

https://epubs.siam.org/doi/10.1137/0714068 offers an A-L stable 2nd order SDIRK method: SDIRK2

-----

>= 3 stage methods : https://www.sciencedirect.com/science/article/pii/S002199912400281X

Start by reading the linked DIRK and SDIRK2 references, then locate MESA's current backward-Euler time-integration and Jacobian evaluation entry points. Determine whether a two-stage DIRK method is numerically appropriate and define convergence and stability tests before implementation; done means the optional method works without regressing the default solver.

Written by the indexing model from the issue text.