This paper presents dynamic modeling of rotor bearing systems with rigid disks, distributed parameter finite rotor elements and flexible, discrete multibearings. The previous works have included the effects of rotary inertia, gyroscopic moment, axial load, internal viscous and hysteretic damping and transverse shear deformations, but have not considered them all in the same model. A computer program is developed in this work to calculate the forward and backward whirl speeds, the corresponding mode shapes, the dynamic unbalance response of multibearing rotor systems and to evaluate rotor stability. It utilizes the banded property of the system matrices to reduce the computational effort for the complex eigensolution. A combined bisection and inverse iteration technique is used for the complex eigen-value problem. The combined effect of shear deformations and internal damping on the stability, forward and backward whirl speeds, as well as on the response to unbalance excitation, are investigated with the numerical examples given. This study confirms that any inherent material damping in the shaft itself does not affect the unbalance response if the bearings are isotropic, which is, however, not the case when the bearings are orthotropic causing elliptical synchronous whirl orbits.