In the paper, we discuss the regular fractional Sturm-Liouville problem in a bounded domain, subjected to the homogeneous mixed boundary conditions. The results on exact and numerical solutions are based on transformation of the differential fractional Sturm-Liouville problem into the integral one. First, we prove the existence of a purely discrete, countable spectrum and the orthogonal system of eigenfunctions by using the tools of Hilbert-Schmidt operators theory. Then, we construct a new variant of the numerical method which produces eigenvalues and approximate eigenfunctions. The convergence of the procedure is controlled by using the experimental rate of convergence approach and the orthogonality of eigenfunctions is preserved at each step of approximation. In the final part, the illustrative examples of calculations and estimation of the experimental rate of convergence are presented.