Should sol.u be a Gaussian? #293
Description
This would make it more clear what the solvers are actually doing. But, it might also make the solvers slightly less compatible with any code that just uses sol.u. I think overall we might want to do this, as then it's more consistent as all the solution values are Gaussians, whereas right now there's sol.u which is just the mean, sol.pu the marginals, but interpolations sol(t) actually return a Gaussian again. If everything is Gaussian, it's more consistent.