Abstract This work presents a function interpolation method based on a Runge–Kutta (RK) scheme to achieve smooth tracking of periodic solutions and their bifurcation behaviors in nonlinear dynamical systems. Unlike conventional harmonic balance (HB) methods that rely on the Newton–Raphson (NR) iteration, the proposed method converts nonlinear equations into linear equations with periodic coefficients, and advances solutions along tangent directions determined by linear equations using RK-based interpolation. This framework enables robust tracking without requiring iterative correction or carefully tuned initial guesses,especially near bifurcation points. A systematic bifurcation analysis is performed by Hill’s method based on the triangular collocation method, allowing explicit identification of tangent directions associated with saddle-node (SN), symmetry-breaking, period-doubling, and Neimark–Sacker (NS) bifurcations. As a result, the method achieves smooth transitions between main branches and bifurcation branches within a unified numerical framework. The effectiveness and generality of the proposed approach are demonstrated through three representative nonlinear systems: the coupled van der Pol equation, the Mathieu–Duffing equation, and the van der Pol–Mathieu equation with external excitation. Numerical results show that the FIM can provide accurate bifurcation tracking with improved robustness near critical points while improves tracking performance and reduces computational cost relative to the incremental HB method by arc-length increments or the arc-length continuation method.
Zhang et al. (Wed,) studied this question.