In this paper, we discuss a class of eigenvalue problems of fractional differential equations of order $\alpha\in (3, 4]$ with variable coefficients. The method of solution is based on utilizing the fractional series solution to find theoretical eigenfunctions. Then, the eigenvalues are determined by applying the associated boundary conditions. A notable result, for certain cases, is that the eigenfunctions are characterized in terms of the Mittag-Leffler or semi Mittag-Leffler functions. The efficiency and accuracy of the present algorithm are demonstrated through several numerical examples.