The boundary element method (BEM) is developed for nonhomogeneous bodies. The static or steadystate response of such bodies leads to boundary value problems for partial differential equations (PDEs) of elliptic type with variable coefficients. The conventional BEM can be employed only if the fundamental solution of the governing equation is known or can be established. This is, however, out of question for differential equations with variable coefficients. The presented method uses simple, known fundamental solutions for homogeneous isotropic bodies to establish the integral equation. An additional field function is introduced, which is determined from a supplementary domain integral boundary equation. The latter is converted to a boundary integral by employing a domain meshless technique based on global approximation by radial basis function series. Then the solution is evaluated from its integral representation based on the known fundamental solution. The presented method maintains the pure boundary character of the BEM, since the discretization into elements and the integrations are limited only to the boundary. Without restricting its generality, the method is illustrated for problems described by second-order differential equations. Therefore, the employed fundamental solution is that of the Laplace equation. Several problems are studied. The obtained numerical results demonstrate the effectiveness and accuracy of the method. A significant advantage of the proposed method is that the same computer program is utilized to obtain numerical results regardless of the specific form of the governing differential operator.