Summary: We investigate the numerical solution of regular and singular Sturm-Liouville problems by means of finite difference schemes of high order. In particular, a set of difference schemes is used to approximate each derivative independently so to obtain an algebraic problem corresponding to the original continuous differential equation. The endpoints are treated depending on their classification and in case of limit points, no boundary condition is required. Several numerical tests are finally reported on equispaced grids show the convergence properties of the proposed approach.