  A matrix method for the solution of direct fractional Sturm-Liouville
problems on bounded domain is proposed where the fractional derivative is
defined in the Riesz sense. The scheme is based on the application of the
Galerkin spectral method of orthogonal polynomials. The order of convergence of
the eigenvalue approximations with respect to the matrix size is studied. Some
numerical examples that confirm the theory and prove the competitiveness of the
approach are finally presented.
