In this article, a new operational matrix of fractional integration of Hermite polynomials is derived to solve multi-order linear fractional differential equations (FDEs) with spectral tau approach. We firstly convert the FDEs into an integrated-form through multiple fractional integration in association with the Riemann-Liouville sense. This integral equation is then for-mulated as an algebraic equation system with Hermite polynomials. Finally, linear multi-order FDEs with initial conditions are solved with this method. We present exact and approximated solutions for a number of representative examples. Numerical results indicate that the pro-posed method provides a high degree of accuracy to solve the linear multi-order FDEs.