The conventional boundary element method (BEM) uses internal cells for the domain integrals, when nonlinear problems or problems with domain effects are solved. In the conventional BEM, however, the merit of the BEM, which is easy preparation of data, is lost. This paper presents numerical integration for a meshless BEM, which does not require the internal cells. This method uses arbitrary internal points instead of internal cells. In this paper, a multi-dimensional interpolation method is proposed. Furthermore, multidimensional numerical integration is proposed using the above multidimensional interpolation method. This method is based on an improved multiple-reciprocity BEM (triple-reciprocity BEM) for heat conduction analysis with heat generation. In order to investigate the efficiency of this method, several numerical examples are given.