A dual-level method of fundamental solutions in conjunction with kernel- independent fast multipole method is proposed in this study. The competitive attributes of the method are that it inherits high accuracy of the method of fundamental solutions, yet avoids producing the resulting ill-conditioned linear system of equations. In contrast to the method of fundamental solutions, the proposed method places two sets of source nodes on the fictitious boundary and physical boundary, respectively, and then combines the fundamental solutions generated by these two sets of source nodes as the modified fundamental solutions of the Laplace equation. This strategy improves significantly the stability of the method of fundamental solutions. In addition, the method is accelerated by the kernel-independent fast multipole method, which reduces the asymptotic complexity of the method to $\mathcal{O}(N)$ from $\mathcal{O}({N}^{2})$. Numerical experiments show that the method can simulate successfully the large-scale heat conduction problems via a single laptop with up to 250000 degrees of freedom.