Implement a preconditioned version as suggested by Krämer et al.

[nathanaelbosch]

The paper "Stable Implementation of Probabilistic ODE Solvers" by Krämer and Hennig (see arXiv) suggests preconditioned state transition matrices which then do not depend on the step-size, together with a diagonal preconditioner that does. Is it possible to compute these quantities directly in square-root form? If yes, that would be a great addition to the package.