This paper is devoted to study a computation scheme to approximate solution of fractional differential equations (FDEs) and coupled system of FDEs with variable coefficients. We study some interesting properties of shifted Legendre polynomials and develop a new operational matrix. The new matrix is used along with some previously derived results to provide a theoretical treatment to approximate the solution of a generalized class of FDEs with variable coefficients. The new method have ability to convert fractional order differential equations having variable coefficients to system of easily solvable algebraic equations. We gave some details to show the convergence of the scheme. The efficiency and applicability of the method is shown by solving some test problems. To show high accuracy of proposed method we compare out results with some other results available in the literature. The proposed method is computer oriented. We use M atLab to carry out necessary calculations.