Buckling analysis of orthotropic graphene monolayer and multilayer nanoplates under bi-axial compression were studied in this work. Graphene has extraordinary properties such as nano-scale size, high mechanical stiffness and roughness, high electrical and thermal conductivity, flexibility (ductility) and magnetic property. The unique properties of this material have received a lot of attention from the scientists. Because of the importance of buckling in these plates, thermal buckling of the graphene sheets has been studied. derivation of constitutive equations, the small scale effects were taken into account using the Eringen’s nonlocal elasticity theory. This theory was employed in this study because the results obtained from this theory have been proved to be in good agreement with those obtained from atomic simulations. The governing equations for the orthotropic nano-plates buckling problem were derived using the classic theory of plates. In order to solve the obtained differential equations, the so-called Galerkin method was employed. Solving eigenvalue problems, like current buckling one, can be considered as one of the strength point of Galerkin method in numerical solution. The solution of algebraic equations in one and also two stages at the one-stage and two-stage analytical format, which made it possible to study different boundary conditions, was presented. The Lagrangian and harmonic shape functions, are assumed in one and two-stage analysis, respectively. Different problems with various boundary conditions were solved using a developed computer code in the MATLAB software environment, and the effect of various parameters on the critical load was studied. As for the verifications purposes, the results obtained from this work were compared with those obtained from Navier’s method and excellent agreement were found between the two sets of results. For the better conclusio a new parameter called the thermal load ratio was defined. Also the critical buckling load in the X direction was normalized. For single-layered nanoplates, thermal effects on the buckling load were studied for the rhombic and trapezoidal geometry and the variations of orthotropic ratio and also different types of supports such as simply and clamped ones was applied in thermal stability analysis and their effects on the thermal load were investigated. In addition, the effect of both small scale coefficient at different modes and the effective geometrical parameters on the thermal buckling load was studied. In a different type of analysis, a particular support was chosen and the range of plate dimensions in which the Eringen’s nonlocal elasticity is valid, was foundThe investigated effect of geometry parameters such as rhomb angle, trapezoid edges ratio, skew angle and rectangle aspect ratio, on the buckling load was plotted in the form of separate diagrams for different scale coefficient. For multilayer nanoplates, the effects of different parameters such as Winkler’s stiffness, Pasternac’s stiffness and interlayer interactions on thermal buckling variations were also investigated. For the multilayer nanoplates, the Winkler-Pasternac’s model was utilized in the simulation of interactions between the interlayer Vanderwalls’s and the elastic substrate. Our results indicate that each of these parameters is directly proportional to the change in buckling load. Finally, the effects of both the number of layers and the mechanical properties of layeres on stability were discussed. This study shows that an increase in the number of layeres considerably increases in stability range. keywords Thermal stability, Single/multi-layered orthotrpic nanoplate, Nonlocal elasticity,Galerkin method