The traditional Method of Fundamental Solutions is revisited, based on a special strategy of defining the external source points. Unlike the classical Method of Fundamental Solutions, the sources are categorized into groups; the density of the spatial distribution of the sources decreases rapidly far from the boundary. On each group, the original problem is discretized using the same set of boundary collocation points. Such groups of sources are constructed in a fully automated way by the quadtree/octtree algorithm. The discretized problems are solved in the sense of least squares. A simple multi-level method is built up, using the (conjugate) gradient iteration as a smoothing procedure. The resulting method significantly reduces the computational complexity. Moreover, the problem of evaluation singular integrals as well as the problem of severely ill-conditioned matrices are avoided. The method is generalized to 3D axisymmetric potential problems as well.