A new meshless local Petrov-Galerkin (MLPG) approach is introduced for solving linear potential problems with high accuracy. The method uses a local symmetric weak form (LSWF) and moving least squares (MLS) approximation. Essential boundary conditions are imposed via a penalty method, eliminating the need for a finite element mesh for interpolation or integration. All integrals are computed over spheres (or circles) centered at each point, ensuring efficiency and accuracy. The method is simple, efficient, and suitable for both linear and nonlinear problems. It avoids the need for domain or boundary elements, making it more flexible than conventional methods like FEM, EFG, and BEM. The MLS approximation is used to interpolate the unknown variable, with coefficients determined by minimizing a weighted discrete L2 norm. The method is applicable to various problems, including elasticity and multi-dimensional boundary value problems. The MLS approximation is defined by a complete monomial basis, with coefficients determined by solving a linear system. The method is well-defined when the matrix A is non-singular, requiring at least m non-zero weight functions for each sample point. The approach is validated through numerical examples, showing high convergence rates and accuracy for Sobolev norms. The method is efficient and attractive for engineering applications.A new meshless local Petrov-Galerkin (MLPG) approach is introduced for solving linear potential problems with high accuracy. The method uses a local symmetric weak form (LSWF) and moving least squares (MLS) approximation. Essential boundary conditions are imposed via a penalty method, eliminating the need for a finite element mesh for interpolation or integration. All integrals are computed over spheres (or circles) centered at each point, ensuring efficiency and accuracy. The method is simple, efficient, and suitable for both linear and nonlinear problems. It avoids the need for domain or boundary elements, making it more flexible than conventional methods like FEM, EFG, and BEM. The MLS approximation is used to interpolate the unknown variable, with coefficients determined by minimizing a weighted discrete L2 norm. The method is applicable to various problems, including elasticity and multi-dimensional boundary value problems. The MLS approximation is defined by a complete monomial basis, with coefficients determined by solving a linear system. The method is well-defined when the matrix A is non-singular, requiring at least m non-zero weight functions for each sample point. The approach is validated through numerical examples, showing high convergence rates and accuracy for Sobolev norms. The method is efficient and attractive for engineering applications.