This study introduces an improved numerical algorithm that is capable of analyzing nonlinear vibrations and bifurcations of general, finite, large-order rotordynamic systems supported on nonlinear bearings. An industrial rotor generally consists of several sections and stages, but numerical shooting/continuation method has been applied to a simple Jeffcott type rotor instead of complex models due to the computational burden of the numerical procedure; it becomes significant when the rotor combined with nonlinear finite bearing models. Here, some mathematical/computational techniques such as a deflation algorithm and the parallel computing are suggested for acceleration along with the conventional treatment of model reduction scheme. An eight-stage compressor rotor supported by two identical five-pad tilting pad journal bearings (TPJB) is selected as a mechanical model to test the numerical incorporation of the algorithms. The rotor beam is modelled with 35 nodes, 140 DOF based on Euler beam theory, and the fluid reaction forces from the two TPJB are calculated using simplex, triangular type finite meshes on the pads. In the numerical procedure, the shooting/continuation combined with the acceleration schemes identifies the solution curves of periodic responses and determines their stability. The orbital motions of coexistent responses are obtained from the solution manifolds.