Dynamic stiffness equations are formulated for variable thickness cylindrical shells, under the assumptions of Donnell, Timoshenko and Flugge theories. Transcendental dynamic stiffness matrices are formed by solving numerically the governing eighth order differential equations using the boundary-value solver COLSYS. Undamped natural frequencies are found using the Wittrick-Williams algorithm. The shell is divided into smaller elements whose clamped end natural frequencies exceed the highest frequency of interest. A parametric study examines the effects of varying the geometry, degree of thickness taper and end conditions on the natural frequencies and mode shapes, providing benchmark results for designers and future researchers.