A Galerkin boundary node method (GBNM), for boundary only analysis of partial differential equations, is discussed in this paper. The GBNM combines an equivalent variational form of a boundary integral equation with the moving least-squares (MLS) approximations for generating the trial and test functions of the variational formulation. In this approach, only a nodal data structure on the boundary of a domain is required, and boundary conditions can be implemented directly and easily despite of the fact that the MLS shape functions lack the delta function property. Formulations of the GBNM using boundary singular integral equations of the second kind for potential problems are developed. The theoretical analysis and numerical results indicate that it is an efficient and accurate numerical method.