This paper is devoted to studying the variational approximation for the higher order regular fractional Sturm-Liouville problems (FSLPs). Using variational principle, we demonstrate that the FSLP has a countable set of eigenvalues and corresponding unique eigenfunctions. Furthermore, we establish two results showing that the eigenfunctions corresponding to distinct eigenvalues are orthogonal, and the smallest (first) eigenvalue is the minimizer of the functional. To validate the theoretical result, we also present a numerical method using polynomials phi j(t)=tj+1(1-t)2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPhi _j(t) = t<^>{j+1}(1-t)<^>2$$\end{document} for j=1,2,3,MIDLINE HORIZONTAL ELLIPSIS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,2,3,\cdots $$\end{document} as a basis function. Further, the Lagrange multiplier method is used to reduce the fractional variational problem into a system of algebraic equations. In order to find the eigenvalues and eigenfunctions, we solve the algebraic system of equations. Further, the analytical convergence and the absolute error of the method are analyzed.