This paper presents a numerical method for solving a class of fractional variational problems (FVPs) with multiple dependent variables, multi order fractional derivatives and a group of boundary conditions. The fractional derivative in the problem is in the Caputo sense. In the presented method, the given optimization problem reduces to a system of algebraic equations using polynomial basis functions. An approximate solution for the FVP is achieved by solving the system. The choice of polynomial basis functions provides the method with such a flexibility that initial and boundary conditions can be easily imposed. We extensively discuss the convergence of the method and finally present illustrative examples to demonstrate validity and applicability of the new technique.