A new meshless method to solve partial differential equations, especially the Laplace equation, is presented. The method is based on a multi-elliptic interpolation invented in a series of previous papers by the author. The idea is to use homogeneous solutions to a higher order elliptic operator, for example \(\Delta(I-\frac{1}c\Delta)\), for both finding a particular and a homogeneous solution to the problem. The solution procedure involves a quadtree based multi-level method to speed up the solution process. Two variants of the method are presented, one including a free parameter \(c\) and a second parameter free called the pure biharmonic boundary interpolation method. For the pure biharmonic boundary interpolation method the author proves an error estimate for the obtained numerical solution. Convergence studies for the Laplace equation on a disc are provided.