The subject matter of this dissertation is the design, analysis and practical implementation of a new numerical method to approximate the eigenvalues and eigenfunctions of regular Sturm–Liouville problems, given in Liouville’s normal form, defined on compact intervals, with self-adjoint separated boundary conditions. These are classical problems in computational mathematics which lie on the interface between numerical analysis and spectral theory, with important applications in physics and chemistry, not least in the approximation of energy levels and wave functions of quantum systems. Because of their great importance, many numerical algorithms have been proposed over the years which span a vast and diverse repertoire of techniques. When compared with previous approaches, the principal advantage of the numerical method proposed in this dissertation is that it is accompanied by error bounds which: (i) hold uniformly over the entire eigenvalue range, and, (ii) can attain arbitrary high-order. This dissertation is composed of two parts, aggregated according to the regularity of the potential function. First, in the main part of this thesis, this work considers the truncation, discretization, practical implementation and MATLAB software, of the new approach for the classical setting with continuous and piecewise analytic potentials (Ramos and Iserles, 2015; Ramos, 2015a,b,c). Later, towards the end, this work touches upon an extension of the new ideas that enabled the truncation of the new approach, but instead for the general setting with absolutely integrable potentials (Ramos, 2014).