Pile buckling is infrequent, but sometimes it can occur in slender piles (i.e., micropiles) driven into soils with soft layers and/or voids. Buckling analysis of piles becomes more complex if the pile is surrounded by multi-layered soil. In this case, the well-known Timoshenko’s solution for pile buckling is of no use because it refers to single-layered soils. A variational approach for buckling analysis of piles in multi-layered soils is herein proposed. The proposed method allows for the estimation of the critical buckling load of piles in any multi-layered soil and for any boundary condition, provided that the distribution of the soil coefficient of the subgrade reaction is available. An eigenvalue-eigenvector problem is defined, where each eigenvector is the set of coefficients of a Fourier series describing the second-order displaced shape of the pile, and the related buckling load is the eigenvalue, thus obtaining the effective buckling load as the minimum eigenvalue. Besides the pile deformed shape, the stiffness distribution in the multi-layered soil is also described through a Fourier series. The Rayleigh–Ritz direct method is used to identify the Fourier development coefficients describing the pile deformation. For validation, buckling analysis results were compared with those obtained from an experimental test and a finite element analysis available in the literature, which confirmed this method’s reliability.