In this work, the least squares weighted residual method is used to solve the two-dimensional (2D) elasticity problem of a rectangular plate of in-plane dimensions 2a 2b subjected to parabolic edge tensile loads applied at the two edges x = a. The problem is expressed using Beltrami–Michell stress formulation. Airy’s stress function method is applied to the stress compatibility equation, and the problem is expressed as a boundary value problem (BVP) represented by a non-homogeneous biharmonic equation. Airy’s stress functions are chosen in terms of one and three unknown parameters and coordinate functions that satisfy both the domain equations and the boundary conditions on the loaded edges. Least squares weighted residual integral formulations of the problems are solved to determine the unknown parameters and thus the Airy stress function. The normal and shear stress fields are determined for the one-parameter and the three-parameter coordinate functions. The solutions for the stress fields are found to satisfy the stress boundary conditions as well as the domain equation. The presented solutions for the Airy stress function and the normal stresses and shear stress fields are identical with solutions obtained by using variational Ritz methods, Bubnov–Galerkin methods and agree with results obtained by Timoshenko and Goodier.