In this paper we present an efficient numerical algorithm for solving linear and nonlinear systems of fractional differential equations (FDEs). The homotopy analysis method is applied to construct the numerical solutions. The proposed algorithm avoids the complexity provided by other numerical approaches. The method is applied to solve three systems of FDEs. Results obtained using the scheme presented here agree well with the analytical solutions and the numerical results presented elsewhere. Moreover, an attempt has been made to address few issues like the effect of varying the auxiliary parameter ℏ, auxiliary function H(t), and the auxiliary linear operator L on the order of local error and convergence of the method. Also, we show that the homotopy perturbation method and Adomain decomposition method are special cases of the homotopy analysis method. © 2009 Elsevier Inc. All rights reserved.