In this paper, the spectral collocation method is investigated for the numerical solution of multi-order Fractional Differential Equations (FDEs). We choose the orthogonal Jacobi polynomials and set of Jacobi Gauss–Lobatto quadrature points as basis functions and grid points respectively. This solution strategy is an application of the matrix-vector-product approach in spectral approximation of FDEs. The fractional derivatives are described in the Caputo type. Numerical solvability and an efficient convergence analysis of the method have also been discussed. Due to the fact that the solutions of fractional differential equations usually have a weak singularity at origin, we use a variable transformation method to change some classes of the original equation into a new equation with a unique smooth solution such that, the spectral collocation method can be applied conveniently. We prove that after this regularization technique, numerical solution of the new equation has exponential rate of convergence. Some standard examples are provided to confirm the reliability of the proposed method.