This paper presents an advanced application of the dynamic stiffness method for the free vibration analysis of moderately thick circular cylindrical shells, based on a generalized Flügge shell theory accounting for shear deformation, rotary inertia and effects of initial stresses. Unlike previous studies, the governing differential equations are solved exactly for each frequency of interest, eliminating the need for numerical approximations in the solution process. An exact dynamic stiffness matrix derived from the strong-form solution is developed for a fully free cylindrical shell element and implemented in a genuine MATLAB code to efficiently compute natural frequencies and mode shapes. The numerical study includes examples featuring stepwise thickness variations, intermediate supports, and initial stresses, providing insights into a wide range of structural applications. The results are validated through comparison with finite element analysis and published data, demonstrating the accuracy, reliability, and computational efficiency of the proposed approach for complex cylindrical shell structures. Additionally, the proposed method addresses limitations of previous studies by capturing all relevant natural frequencies. Finally, numerous high-accuracy results are provided to serve as benchmark solutions for validating future research in this field.