Complex multi-body systems are widely used in various fields such as transportation, mechanical, aerospace, civil and weapon engineering. This paper proposes an exact, efficient, reliable and general dynamic stiffness (DS) method for analytical modelling multibody systems, where any number of rigid bodies can be connected to any number of flexible beams at any points with any angles. Different DS beam elements are developed based on arrange of different theories, such as the classical, Rayleigh-Love and Mindlin-Herrmann theories for axial vibration, while the Bernoulli-Euler and Timoshenko theories for bending vibration. Therefore, the dynamics of beams with both lower and higher aspect ratios are described exactly within the whole frequency range. Moreover, both the planar and spatial assembly procedures for connecting the beam and rigid bodies with any kind of eccentricity and at any angles is simplified without resorting to matrix inversion. Then, the global dynamic stiffness matrix of multi-body systems is obtained. The difficulty generally encountered in computing the problematic J0 count when applying the Wittrick-Williams algorithm for modal analysis has been overcome. Applications of different beam theories for both axial and bending vibrations have enabled the examination of the role played by rigid-body parameters on the multi-body system’s dynamic behaviour. Finally, some exact benchmark results are provided and compared with published results and with finite element solutions. This research provides an exact and highly efficient analysis tool for multi-body system dynamics which is for the free vibration analysis, ideally suited for optimization and inverse problems and can be extended to the analytical model of beam-plate element coupled rigid bodies structures.