In this paper an energy method is presented for the linear buckling analysis of first order shear deformable plates. The displacement fields are defined in terms of the shape functions, which correspond to a set of predefined points and are composed of significantly high order polynomials. The locations of these points are found by mapping the geometry using the naturalized coordinates and bilinear shape functions. In order to evaluate the method, fully clamped and simply supported rectangular plates subjected to uniform uniaxial compressive loading on two opposite edges of the plate are investigated thoroughly and the results are compared with the exact solution given in the monograph of Timoshenko and Gere (1961). The method is extended to the analysis of perforated plates, wherein the negative stiffness computed over the opening area from in-plane and out-of-plane deformation modes is superimposed to the stiffness of the full plate. Numerical results are then favorably compared with those obtained by finite element methods. Other cases such as; rectangular plates with eccentrically located openings of different shapes are studied and reported in this paper with regards to the effect of aspect ratio, hole size, and hole position on the buckling. For a square plate with a large circular opening at the center, diameter being 80 percent of the length, the present method yields buckling coefficient 12.5 percent higher than the one from the FEM.