In this work, on the basis of homotopy analysis method, we propose a numerical procedure for solving multi-order fractional differential equations by transforming the original multi-order fractional differential equation into a system of single-order equations. To our knowledge this paper presents the first viable numerical method for an analytic solution of multi-order fractional differential equations. This approach provides a simple way to ensure the convergence of series of solution so that one can always get accurate enough approximations.