The Boundary Integral Equation Method reduces the spatial dimension of an elliptic problem by converting the original n-dimensional partial differential equation to an (n - 1)-dimensional boundary integral equation defined on the boundary of the domain of the original problem. At the same time, the discretisation of the problem is also remarkably simplified. The price of these advantages, however, is that the structure, as well as the algebraic properties of the resulting boundary element matrices, are somewhat unpleasant, since they are neither sparse nor self-adjoint in general, even if the original problem is self-adjoint. Consequently, the computational cost of the Boundary Integral Equation Method seems to be unnecessarily high. To make the method more economic from a computational point of view, we present a numerical technique based on the multipole expansion method, which reduces the computational cost of the appearing matrix-vector multiplications (i.e. the evaluations of the discretized boundary integrals) by a remarkable amount. The method is applicable also to die reconstruction problem, when the inner solution is to be reconstructed from the boundary solutions. © 1998 Elsevier Science S.A.