<p>Finite Element Method (FEM) simulations involving elasto-plasticity and inhomogeneous response can be accelerated by the use of biased mesh with higher-order approximation. In this regard, polygonal elements with higher-order shape functions can be useful. This is addressed in the present work where quadratic Wachspress shape functions (quotient of two) for strictly convex polygonal patches are developed. The shape functions are subsequently used in Galerkin Polygonal Finite Element Method (P-FEM) to solve elasto-plastic boundary value problems. The accuracy of the quadratic P-FEM is analyzed by comparing its response with linear P-FEM, analytical solution and conventional FEM for three different problems involving a rectangular domain, a cantilever beam and a square plate with a hole having elastic and elasto-plastic constitutive behaviour. The comparisons show that quadratic P-FEM can capture the local responses accurately with a significantly coarser mesh. The h-refinement of elements based on the reduction of center-to-center distances of polygons also shows an improvement in accuracy, emphasizing the convergent nature of the method.</p>