Homotopy analysis method (HAM) is a popular semianalytical method used widely in applied sciences. It stands out from the rest of the semi-analytical methods as it provides a family of solutions to nonlinear equations, including ordinary differential equations (ODEs), partial differential equations, etc. The convergence characteristics of the solutions can be varied by changing an auxiliary parameter (h) in HAM. The convergence region of solution of ODEs using HAM can be improved by applying it over multiple intervals of time, which is referred to as multi-stage HAM (MHAM). In this paper, MHAM models for the IEEE Model 2.2 synchronous machine, IEEE Type-1 excitation system, first-order governor and first-order turbine models have been developed. The applicability of MHAM for power system dynamic simulations has been investigated in this paper using seven widely used test systems ranging from 10 generators 39 bus systems to 4092 generators 13 659 bus systems. The effect of number of terms, h and the time step on the accuracy and stability of the solution has been studied. The effectiveness of MHAM has been compared with the modified Euler (ME) and midpoint Trapezoidal (TrapZ) methods. The accuracy of MHAM has been found to be comparable with ME and TrapZ methods for the values of h between -1.05 and -0.95. The best accuracy is obtained for h = -1.0, which is a special case of MHAM called multi-stage homotopy perturbation method (MHPM). In this paper, it is also shown that MHPM is equivalent to multi-stage adomian decomposition method, which has been recently explored for large power system simulations.