An efficient domain decomposition method is proposed for solving the free, harmonic and transient vibrations of isotropic and composite cylindrical shells subjected to various combinations of classical and non-classical boundary conditions. Multi-segment partitioning strategy is adopted to accommodate the computing requirements of high-order vibration modes and responses. The continuity constraints on the segment interfaces are incorporated into the system potential functional by means of a modified variational principle and least-squares weighted residual method. An arbitrarily laminated version of Reissner–Naghdi’s shell theory is employed to formulate the theoretical model. Double mixed series, i.e., the Fourier series and orthogonal polynomials, are used as admissible displacement functions for each shell segment. The utility and robustness of the method for the application of various basis functions are evaluated with the following four sets of orthogonal polynomial series, i.e., the Chebyshev orthogonal polynomials of first and second kind, Legendre orthogonal polynomials of first kind and Hermite orthogonal polynomials. To test the convergence, efficiency and accuracy of the present method, free and forced vibrations (including the harmonic and transient vibrations) of isotropic and composite laminated cylindrical shells are examined under different combinations of free, shear-diaphragm, simply-supported, clamped and elastic supported boundaries. The theoretical results are compared with those previously published in literature, and the ones obtained by using the finite element program ANSYS. Very good agreement is observed.