  We propose a fast algorithm for the probabilistic solution of boundary value
problems (BVPs), which are ordinary differential equations subject to boundary
conditions. In contrast to previous work, we introduce a Gauss--Markov prior
and tailor it specifically to BVPs, which allows computing a posterior
distribution over the solution in linear time, at a quality and cost comparable
to that of well-established, non-probabilistic methods. Our model further
delivers uncertainty quantification, mesh refinement, and hyperparameter
adaptation. We demonstrate how these practical considerations positively impact
the efficiency of the scheme. Altogether, this results in a practically usable
probabilistic BVP solver that is (in contrast to non-probabilistic algorithms)
natively compatible with other parts of the statistical modelling tool-chain.
