This study investigated the elastic static stability analysis of homogeneous and isotropic thick rectangular plates with twelve boundary conditions and carrying uniformly distributed uniaxial compressive load using the direct variational method. In the analysis, a thick plate energy expression was developed from the three-dimensional (3-D) constitutive relations and kinematic deformation; thereafter the compatibility equations used to resolve the rotations and deflection relationship were obtained. Likewise, the governing equations were derived by minimizing the equation for the potential energy with respect to deflection. The governing equation is solved to obtain an exact deflection function which is produced by the trigonometric and polynomial displacement shape function. The degree of rotation was obtained from the equation of compatibility which when equated to the deflection function and put into the potential energy equation formulas for the analysis were obtained after differentiating the outcome with respect to the deflection coefficients. The result obtained shows that the non-dimensional values of critical buckling load decrease as the length-width ratio increases (square plate being the highest value), this continues until failure occurs. This implies that an increase in plate width increases the probability of failure in a plate. Hence, it can be deduced that as the in-plane load on the plate increase and approaches the critical buckling, the failure in a plate structure is abound to occur. Meanwhile, the values of critical buckling load increase as the span-thickness ratio increases for all aspect ratios. This means that, as the span-thickness ratio increases an increase in the thickness increases the safety in the plate. It also indicates that the capacity of the plate to resist buckling decreases as the span-depth ratio increases. To establish the credibility of the present study, classical plate theory (CPT), refined plate theory (RPT) and exact solution models from different studies were employed to validate the results. The present works critical buckling load varied with those of CPT and RPT with 7.70% signifying the coarseness of the classical and refined plate theories. This amount of difference cannot be overlooked. The average total percentage differences between the exact 3-D study (Moslemi et al., 2016), and the present model using polynomial and trigonometric displacement functions is less than 1.0%. These differences being so small and negligible indicates that the present model using trigonometric and polynomial produces an exact solution. Thus, confirming the efficacy and reliability of the model for the 3-D stability analysis of rectangular plates.