It is very difficult to obtain perfect analytical solutions of differential equations of gear system dynamics. The dynamic model which describes the torsional vibration behaviors of gear system has been introduced accurately in this paper. The differential equation of gear system nonlinear dynamics exhibiting combined nonlinearity influence such as time-varying stifffiess, tooth backlash and dynamic transmission error (DTE) has been proposed by using a polynomial of degree 7 to fit the nonlinear backlash function and using time-varying stiffness Fourier transform to obtain its harmonic forms with 5 orders for the first time. The theory of GEAR method has been presented. Contrasting with other numerical methods, GEAR method has higher precision and calculation efficiency, especially in solving stiff differential equations. Based on GEAR method, the numerical calculation method for solving differential equations of gear system dynamics has been presented in this paper. Numerical calculation results proved that the numerical solutions by using Gear method is validated by comparison with experimental measurements and can be used to solve all kinds of differential equations, especially for large differential equations.