The approximation properties of shifted fundamental solutions are stated for a general class of second order linear elliptic differential operators. Motivated by the method of the fundamental solutions, the singularities of the fundamental solutions are located in the vicinity of a fixed surface, which is the boundary of a 3-dimensional Lipschitz domain. The approximation is estimated first in Ck-norm. The corresponding convergence rate is determined depending on the density of the singularities. This delivers a solid theoretical basis for the method of fundamental solutions. The results are extended to exterior domains and as a new application, the numerical solution of Dirichlet-to-Neumann problems are constructed and analyzed on unbounded domains. The efficiency of this new method is demonstrated in series of numerical experiments.