A special regularization method based on higher-order partial differential equations is presented. Instead of using the fundamental solution of the original partial differential operator with source points located outside of the domain, the original second-order partial differential equation is approximated by a higher-order one, the fundamental solution of which is continuous at the origin. This allows the use of the method of fundamental solutions (MFS) for the approximate problem. Due to the continuity of the modified operator, the source points and the boundary collocation points are allowed to coincide, which makes the solution process simpler. This regularization technique is generalized to various problems and combined with the extremely efficient quadtreebased multigrid methods. Approximation theorems and numerical experiences are also presented.