Multiple solution problems for nonlinear rotor bearing systems (RBSs) are of great interest in research and industrial development; it is a challenge to detect multiple solution distributions for RBSs, especially for large-scale systems. This paper develops a heuristic search method for detecting multiple periodic solution branches of RBSs. In this method, differential equations of RBSs are first converted to algebraic equations by the harmonic balance (HB) method; a model reduction technique is employed to keep only nonlinear degrees of freedom (DOFs). A heuristic search method is developed to detect multiple solutions for RBSs, and multiple solution branches are obtained by numerical continuation. By applying the methodology to a four-DOF lumped parameter system and a rolling bearing rotor system, it is shown that the proposed numerical method can be employed to search multiple solution branches of RBSs efficiently. Meanwhile, isolated solution branch phenomena are observed.