Abstract: The perturbation methods for nonlinear vibration systems make it necessary to solve a set of second-order ordinary differential equations (ODEs), which are obtained by equating the like power of the perturbation solutions respectively. One of the main drawbacks of the ODEs-based methods is of low efficiency, especially for nonlinear vibration systems of multiple degrees of freedom. In this paper, a method of polynomial vectors for solving the approximate solution of nonlinear vibration systems is proposed. The second order ordinary differential equations are written in a set of state equations of the first order first, wherein the nonlinear terms of the state equations are expressed as the products of a constant matrix and a polynomial vector with the like power. By using the direct perturbation method, the linear non-homogeneous equations are obtained for the like power approximations, while the nonlinear terms are written as the products of the constant matrix and the polynomial vector with the previous approximate solutions as its element. Furthermore, the multiplication of polynomials in the polynomial vector is expressed in matrix form via Toeplitz matrix, and then all approximate analytical formulas of the state equations are determined by the first-order non-homogeneous equations. Results show that the proposed method based upon the state equations yields a concise calculation for nonlinear vibration systems of multiply degrees of freedom.