We consider a multi–level method of fundamental solutions for solving polyharmonic problems governed by ΔNu=0,N∈N∖{1} in both two and three dimensions. Instead of approximating the solution with linear combinations of N fundamental solutions, we show that, with appropriate deployments of the source points, it is possible to employ an approximation involving only the fundamental solution of the operator ΔN. To determine the optimal position of the source points, we apply the recently developed effective condition number method. In addition, we show that when the proposed technique is applied to boundary value problems in circular or axisymmetric domains, with appropriate distributions of boundary and source points, it lends itself to the application of matrix decomposition algorithms. The results of several numerical tests are presented and analysed.