<p>This paper presents a combined numerical iteration scheme for efficiently auto-tracking the frequency-response curves of the high-dimensional nonlinear systems, on the basis of the combination of harmonic balance method, pseudo arclength method, and Newton/Lagrange extrapolation techniques. The auto-tracking tangent vectors and auto-adaptive arclengths are introduced. This combined iteration scheme is extended to the complex nonlinear systems with nonlinear damping and time delay; the generalized increment matrix equations and universal flow chart to establish iteration scheme are given. The relationship between the incremental harmonic balance method and harmonic balance-Galerkin-Newton Raphson method is discussed through comparing the separated procedures. On the basis of the harmonic balance-based solutions, the methods to analyze the stability of the steady state solutions and compute the bifurcation points for the complex nonlinear system are presented as well. The combined iteration scheme and the flow chart are applied to the nonlinear differential equations of the systems with high dimensions, time delay and nonlinear damping, respectively. The accuracy of the proposed iteration scheme is demonstrated through its validation with respect to Runge-Kutta solutions. The high computational efficiency of the combined iteration scheme is proved via numerical examples, especially for the high dimensional nonlinear systems.</p>