can we operator split within VODE?

[zingale]

Right now we solve the system:

$X^{n+1}_k = \dot{\omega}_k (\rho, T^{n+1}, X_k^{n+1})$
$e^{n+1} = \dot{S} (\rho, T^{n+1}, X_k^{n+1})$

but what if we treated the energy/temperature update explicitly. We are substepping in VODE, so we could still be accurate / stable (?) if we did:

$X^{n+1}_k = \dot{\omega}_k (\rho, T^{n}, X_k^{n+1})$
$e^{n+1} = \dot{S} (\rho, T^{n}, X_k^{n+1})$

the advantage of this is that the Jacobian then does not need any temperature derivatives, significantly simplifying it and making issues like #729 a non-issue.

[zingale]

we could even do a second order update for e -- in some sense, this becomes an IMEX scheme?