A modified two-timescale incremental harmonic balance (IHB) method is introduced to obtain quasi-periodic responses of nonlinear dynamic systems with combinations of two incommensurate base frequencies. Truncated Fourier coefficients of residual vectors of nonlinear algebraic equations are obtained by a frequency mapping-fast Fourier transform procedure, and complex two-dimensional (2D) integration is avoided. Jacobian matrices are approximated by Broyden's method and resulting nonlinear algebraic equations are solved. These two modifications lead to a significant reduction of calculation time. To automatically calculate amplitude–frequency response surfaces of quasi-periodic responses and avoid nonconvergent points at peaks, an incremental arc-length method for one timescale is extended for quasi-periodic responses with two timescales. Two examples, Duffing equation and van der Pol equation with quadratic and cubic nonlinear terms, both with two external excitations, are simulated. Results from the modified two-timescale IHB method are in excellent agreement with those from Runge–Kutta method. The total calculation time of the modified two-timescale IHB method can be more than two orders of magnitude less than that of the original quasi-periodic IHB method when complex nonlinearities exist and high-order harmonic terms are considered.