The analysis of whole engine rotordynamic models is an important element in the design of aerojet engines. The models include gyroscopic effects and allow for rubbing contact between rotor and stator components such as bladed discs and casing. Due to the nonlinearities inherent to the system, bifurcations in the frequency response may arise. Reliable and efficient methods to determine the bifurcation points and solution branches are required. For this purpose, a multi-harmonic balance approach is presented that allows a numerically efficient detection of bifurcation points and the calculation of both continuous and isolated branches of the frequency response functions. The method is applied to a test case derived from a commercial aeroengine. A bifurcation structure with continuous and isolated solution branches is observed and studied in this paper. The comparison with time marching based on simulations shows both accuracy and numerical efficiency of the newly developed approach.