We consider boundary value problems of the formwith f , g and h continuous in [0, 1]. We use the Liouville-Neumann expansion to design a succession of functions y n (x) that converge uniformly on [0, 1] to the solution y(x) of that problem. In particular, when f (x), g(x) and h(x) are polynomials, y n (x) are also polynomials. We show that the Liouville-Neumann algorithm may be used to approximate eigenvalues and eigenfunctions of eigenvalue problems of the form y + f (x)y + [g(x) + λσ (x)]y = 0, x ∈ (0, 1), y(0) = y(1) = 0, with f , g and σ continuous in [0, 1]. It may be also used to approximate zeros of solutions of regular second-order linear differential equations and, in particular, of some special functions.