This paper analytically studies the buckling of a cylindrical shell having varying thickness under non-uniform axial compressive loads for the first time, which widely exists in engineering practice. A novel quadratic perturbation technique is developed to establish general buckling load formulas for the shell. This method overcomes the difficulties of traditional energy methods in solving high order determinants and deriving direct expressions for buckling loads when shell thickness and axial load are unknown. Applying presented formulas, various shell thicknesses and axial loads are analyzed, and a series of new results for buckling loads are obtained and validated. Even for classical cosine thickness variation under uniform axial compression, we also give general conclusions compared with Koiter’s results by the energy method. The effects of thickness variations and load distribution parameters on buckling loads are analyzed in detail. The presented study in this paper fills the gap and establishes a foundation of buckling analysis for non-uniformly loading cylindrical shells with variable thickness. Certainly, the established formulas are general and available for buckling resistance capacity evaluation for the shells under all circumstances involving thickness variations or/and non-uniform axial compressive loads.