Massive vacuum integrals with a single mass scale are a class of Feynman integrals that appear in many precision calculations within the Standard Model of particle physics and have been calculated to the 4-loop level. In this thesis I start pushing this limit to 5 loops by considering the subclass of fully massive vacuum integrals, which can be used to determine the β-function of Quantum Chromodynamics (QCD). To this end I employ a method devised by Laporta for the evaluation of multi-loop Feynman integrals based on difference equations and factorial series. Significant improvements to this method are introduced to account for the great increase in complexity when going from 4 to 5 loops. An implementation of the improved approach in C ++ is then used to obtain high-precision numerical results for the integrals needed for the 5-loop correction to the QCD β-function.