In this paper, we construct numerical schemes based on the Lagrange polynomial interpolation to solve Fractional Sturm–Liouville problems (FSLPs) in which the fractional derivatives are considered in the Caputo sense. First, we convert the differential equation with boundary conditions into integral form and discretize the fractional integral to generate a system of algebraic equations in the matrix form. Next, we calculate the set of approximate eigenvalues and corresponding eigenvectors. The eigenfunctions are approximated and some of their properties are investigated. The experimental rate of convergence of numerical calculations for the eigenvalues is reported and the order convergence of the numerical method is obtained. Finally, some examples are presented to illustrate the efficiency and accuracy of the numerical method.