The two-dimensional scattered data interpolation problem is investigated. In contrast to the traditional Method of Radial Basis Functions, the interpolation problem is converted to a higher order (biharmonic or modified bi-Helmholtz) partial differential equation supplied with usual boundary conditions as well as pointwise interpolation conditions. To solve this fourth-order problem, the Method of Fundamental Solutions is used. The source points, which are needed in the method, are located partly in the exterior of the domain of the corresponding partial differential equation and partly in the interpolation points. This results in a linear system with possibly large and fully populated matrix. To make the computations more efficient, a localization technique is applied, which splits the original problem into a sequence of local problems. The system of local equations is solved in an iterative way, which mimics the classical overlapping Schwarz method. Thus, the problem of large and ill-conditioned matrices is completely avoided. The method is illustrated via a numerical example.