  Searching for numerical methods that combine facility and efficiency, while
remaining accurate and versatile, is critical. Often, irregular geometries
challenge traditional methods that rely on structured or body-fitted meshes.
Meshless methods mitigate these issues but oftentimes require the weak
formulation which involves defining quadrature rules over potentially intricate
geometries. To overcome these challenges, we propose the Least Squares
Discretization (LSQD) method. This novel approach simplifies the application of
meshless methods by eliminating the need for a weak formulation and
necessitates minimal numerical analysis. It offers significant advantages in
terms of ease of implementation and adaptability to complex geometries. In this
paper, we demonstrate the efficacy of the LSQD method in solving elliptic
partial differential equations for a variety of boundary conditions,
geometries, and data layouts. We monitor h-P convergence across these
parameters and construct an a posteriori built-in error estimator to establish
our method as a robust and accessible numerical alternative.
