  The Method of Fundamental Solutions (MFS) is a popular tool to solve Laplace
and Helmholtz boundary value problems. Its main drawback is that it often leads
to ill-conditioned systems of equations. In this paper we investigate for the
interior Helmholtz problem on analytic domains how the singularities (charge
points) of the MFS basis functions have to be chosen such that approximate
solutions can be represented by the MFS basis in a numerically stable way. For
Helmholtz problems on the unit disc we give a full analysis which includes the
high frequency (short wavelength) limit. For more difficult and nonconvex
domains such as crescents we demonstrate how the right choice of charge points
is connected to how far into the complex plane the solution of the boundary
value problem can be analytically continued, which in turn depends on both
domain shape and boundary data. Using this we develop a recipe for locating
charge points which allows us to reach error norms of typically 10^{-11} on a
wide variety of analytic domains. At high frequencies of order only 3 points
per wavelength are needed, which compares very favorably to boundary integral
methods.
