The finite element analysis for rotor bearing systems has been an essential tool for design, identification, and diagnosis of rotating machinery. Among others, the unbalance response analysis is fundamental in the vibration analysis of rotor bearing systems because rotating unbalance is recognized as a common sourve of vibration in rotating machinery. However there still remains a problem in the aspect of computational efficiency for unbalance response analysis of large rotor bearing systems. Gyroscopic terms and local bearing parameters in rotor bearing systems often make matters worse in unbalance response computation due to the complicated dynamic properties such as rotational speed dependency and/or anisotropy. The present paper proposes an efficient method for unbalance responses of multi-span rotor bearing systems. An improved substructure synthesis scheme is introduced which makes it possible to compute unbalance responses of the system by coupling unbalance responses of substructures that are of self adjoint problem with small order matrices. The present paper also suggests a scheme to easily deal with gyroscopic tems and local, coupling or bearing parameters. The proposed method causes no errors even though the computational effort is reduced drastically. The present method is demonstrated through three test examples.