Although rotors are simplified to be axisymmetric in rotordynamic models, many rotors in the industry are actually non-axisymmetric. Several authors have proposed methods using 3D finite element, rotordynamic models, but more efficient approaches for handling a large number of degrees-of-freedom (DOF) are needed. This task becomes particularly acute when considering parametric excitation that results from asymmetry in the rotating frame. This paper presents an efficient rotordynamic stability approach for non-axisymmetric rotor-bearing systems with complex shapes using three-dimensional solid finite elements. The 10-node quadratic tetrahedron element is used for the finite element formulation of the rotor. A rotor-bearing system, matrix differential equation is derived in the rotor-fixed coordinate system. The system matrices are reduced by using Guyan reduction. The current study utilizes the Floquet theory to determine the stability of solutions for parametrically excited rotor-bearing systems. Computational efficiency is improved by discretization and parallelization, taking advantage of the discretized monodromy matrix of Hsu's method. The method is verified by an analytical model with the Routh–Hurwitz stability criteria, and by direct time-transient, numerical integration for large order models. The proposed and Hill's methods are compared with respect to accuracy and computational efficiency, and the results indicate the limitations of Hill's method when applied to 3D solid rotor-bearing systems. A parametric investigation is performed for an asymmetric Root's blower type shaft, varying bearing asymmetry and bearing damping.