Evaluation of addition coefficients introduced by the addition theorems for vector spherical harmonics is one of the most intractable problems in electromagnetic scattering by multi-sphere systems. The derivation of the analytical expressions for the addition coefficients is lengthy and complex while the computation of the addition coefficients is annoyingly time-consuming even with the reasonably fast computers available nowadays. This paper presents an efficient algorithm for calculating addition coefficients which is based on the recursive relations of scalar addition coefficients. Numerical results from the formulation derived in this paper agree with those of previous published results but the algorithm proposed here reduces the computational time considerably. This paper also discusses the strengths and limitations of other formulations and numerical techniques found in the literature.