Buckling of multilayer graphene sheets (MLGSs) subjected to an axial compressive load in plane-strain condition is studied. Closed-form solutions for buckling load of MLGSs are obtained based on a continuum model for MLGSs. Two different kinematic assumptions, which lead to MLGS beam, which was recently proposed by the authors, and the Euler beam, are used to obtain the buckling loads. The obtained solutions yield significantly different buckling loads when the axial length is small. To validate obtained results, molecular dynamics (MD) simulations are conducted, and they show that the MLGS beam model well captures the buckling load of MLGSs. The buckling solution of MLGS beam model provides two interesting facts. First, the buckling load of MLGSs coincides with the Euler buckling load when the length is large. Second, when the number of layers is large, the buckling strain converges to a finite value, and could be expressed as a linear combination of the buckling strain of single-layer graphene and the ratio between the shear rigidity of interlayer and the tensile rigidity of graphene layer. We validate the asymptotic behavior of buckling strain through MD simulations and show that buckling occurs even when the overall thickness is larger than the axial length. Finally, we present a diagram that contains buckling strain of MLGSs according to the boundary conditions, the number of layers, and the axial length.