The drive axle is the core of the powertrain in automotive vehicles, and numerical simulation is a common and effective approach for designing and analysing drive axles. This article develops a new numerical simulation method for drive axles considering the nonlinear couplings between gears and bearings. Two types of models are proposed and integrated: a reduced model with nonlinear bearing elements and equivalent gear elements, which is efficient for nonlinear bearing simulation, and a gear contact analysis model with equivalent bearing stiffnesses, which can accurately reflect the gear mesh conditions and system deflection. The equivalent gear mesh parameters used in the reduced model are derived from the contact analysis model, and the linear equivalent bearing stiffness matrices used in the contact analysis model are derived from the reduced model. Therefore, the static equilibrium condition of the entire drive axle system is obtained through an iterative calculation process involving the two models under the same load condition. The good correlation between the numerical and experimental results of the gear-loaded contact pattern indicates that the proposed method is reliable.