Krylov methods in truncated EK1

[pnkraemer]

The complexity-dominant operation in the truncated EK1 is a calculation

jax.scipy.linalg.solve_cholesky((S_sqrtm, True), z)

which is used to update the mean.

Sparse Jacobians will imply sparse S_sqrtm (or at least sparse S), and a Krylov method like CG might make the whole shebang faster.

[pnkraemer]

Radau implementation in scipy switches between LU and sparse LU. Not sure whether we wanna do the same, but it might be a useful reference.
https://github.com/scipy/scipy/blob/v1.7.1/scipy/integrate/_ivp/radau.py#L179-L535