  We describe a method to numerically compute multi-loop integrals, depending
on one dimensionless parameter $x$ and the dimension $d$, in the whole
kinematic range of $x$. The method is based on differential equations, which,
however, do not require any special form, and series expansions around singular
and regular points. This method provides results well suited for fast numerical
evaluation and sufficiently precise for phenomenological applications. We apply
the approach to four-loop on-shell integrals and compute the coefficient
function of eight colour structures in the relation between the mass of a heavy
quark defined in the $\overline{\rm MS}$ and the on-shell scheme allowing for a
second non-zero quark mass. We also obtain analytic results for these eight
coefficient functions in terms of harmonic polylogarithms and iterated
integrals. This allows for a validation of the numerical accuracy.
