A new variational formulation for boundary node method (BNM) using a hybrid displacement functional is presented here. The formulation is expressed in terms of domain and boundary variables, and the domain variables are interpolated by classical fundamental solution; while the boundary variables are interpolated by moving least squares (MLS). The main idea is to retain the dimensionality advantages of the BNM, and get a truly meshless method, which does not require a 'boundary element mesh', either for the purpose of interpolation of the solution variables, or for the integration of the 'energy'. All integrals can be easily evaluated over regular shaped domains (in general, semi-sphere in the 3-Dproblem) and their boundaries. Numerical examples presented in this paper for the solution of Laplace's equation in 2-D show that high rates of convergence with mesh refinement are achievable, and the computational results for unknown variable are most accurate. No further integrations are required to compute the unknown variables inside the domain as in the conventional BEM and BNM. Copyright