In this paper, a new reliable algorithm based on an adaptation of the standard homotopy analysis method (HAM) is presented, which is the multistage homotopy analysis method (MSHAM). The freedom of choosing the auxiliary linear operator and the auxiliary parameter are still present in the MSHAM. The solutions of the non-chaotic and the chaotic Chen system which is a three-dimensional system of ordinary differential equations with quadratic nonlinearities were obtained by MSHAM. Numerical comparisons between the MSHAM and the classical fourth-order Runge–Kutta (RK4) numerical solutions reveal that the new technique is a promising tool for solving the non-linear chaotic and non-chaotic Chen system.