We herein propose an algorithm for tracking smooth bifurcations of nonlinear systems with interdependent parameters. The approach is based on a complex formulation of the well-known Harmonic Balance Method (HBM). Hill's method is used to assess the stability of the computed forced response curves and a minimally extended system is built to allow for the parametric continuation of the detected bifurcation points. The feasibility of coupling HBM-based minimally extended systems and arclength continuation algorithms is established and demonstrated. The method offers an efficient way of determining the stability regions of the system. The methodology is applied on a spur gear pair model including the backlash nonlinearity and subjected to transmission error and mesh stiffness fluctuation whose harmonic contents depend on several parameters that do not appear explicitly in the equations of motion.