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 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 non-convergent points at peaks, an incremental arc-length method for one time-scale is extended for quasi-periodic responses with two timescales. Two examples: a Duffing equation and a 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.