The direct application of the classical method of fundamental solutions (MFS) is restricted to homogeneous linear partial differential equations (PDEs). The use of fundamental solutions with different frequencies allowed the extension of the MFS to non-homogeneous PDEs, in particular, for Poisson or Helmholtz equations and for elastostatic or elastodynamic problems. This method has been called method of fundamental solutions for domains (MFS-D), but it faces an approximation problem when the non-homogeneous term presents discontinuities, because the fundamental solutions are analytic functions outside the source point set. In this paper we analyze two domain decomposition techniques for overcoming this approximation problem. The problem is set in the context of the modified Helmholtz equation, and we also establish the missing density results that justify both the MFS and the MFS-D approximations. Numerical results are presented comparing a direct and an iterative domain decomposition technique, with simulations in non-trivial domains.