An exact Fourier series method is developed for the vibration analysis of multispan beam systems. In this method, the displacement on each beam is expressed as a Fourier series expansion plus an auxiliary closed-form function such as polynomials. The auxiliary function is used to deal with all the possible discontinuities, at the end points, with the original displacement function and its derivatives when they are periodically extended over the entire -axis as implied by a Fourier series representation. As a result, not only is it always possible to expand the beam displacements into Fourier series under any boundary conditions, but also the series solution will be substantially improved in terms of its accuracy and convergence. Mathematically, the current Fourier series expansion represents an exact solution to a class of beam problems in the sense that both the governing equations and the boundary/coupling conditions are simultaneously satisfied to any specified degree of accuracy. In the multispan beam system model, any two adjacent beams are generally connected together via a pair of linear and rotational springs, allowing a better modeling of many real-world joints. Each beam in the system can also be independently and elastically restrained at its ends so that all boundary conditions including the classical homogeneous boundary conditions at the end and intermediate supports can be universally dealt with by simply varying the stiffnesses of the restraining springs accordingly, which does not involve any modification of basis functions, formulations, or solution procedures. The excellent accuracy and convergence of this series solution is demonstrated through numerical examples.