A boundary interpolation technique is introduced based on multi-elliptic partial differential equations. The interpolation problem is converted to a special higher order partial differential equation which is completely independent of the geometry of the original problem. Based on this interpolation method, meshless methods are constructed for the 2D Laplace-Poisson equation. The presented approach makes it possible to avoid solving large and dense interpolation equations. The auxiliary higher order partial differential equation is solved by robust, quadtree-based multi-level methods. The results can be easily generalized to 3D problems as well. © 2004 Elsevier Ltd. All rights reserved.