<p>The incremental harmonic balance (IHB) method is commonly used to obtain periodic solutions of strongly nonlinear systems and to track system responses. However, the Newton–Raphson iteration within the IHB method places high demands on initial values, and convergence may not be achieved when these values are not suitable. While initial values can be refined using methods such as the perturbation method and the fast Fourier transform method for a single solution, these approaches are computationally intensive and may not be practical for response tracking. To address convergence issues during response tracking, the Runge–Kutta (RK) method is used here to replace the iteration. In this approach, the IHB method is first applied to obtain a solution at an arbitrary point as the initial value of the RK method. The RK method is then used to track the complete system response, which is referred to as the Runge–Kutta tracking (RKT) method. Using this tracking method, amplitude-parameter and frequency-parameter response curves, as well as phase diagrams, are obtained for the van der Pol oscillator, the coupled van der Pol oscillator, and the oscillator of a slider pendulum under external harmonic excitation. Results are then compared to those from the IHB methods by amplitude increments and frequency increments and the arc-length continuation method with different steps, while stability of the curves is analyzed via the Floquet theory. The system responses from the RKT method agree well with those obtained from the IHB and RK methods, and the RKT method demonstrates better tracking performance.</p>