  The method of the fundamental solutions (MFS) is used to construct an
approximate solution for a partial differential equation in a bounded domain.
It is demonstrated by combining the fundamental solutions shifted to the points
outside the domain and determining the coefficients of the linear sum to
satisfy the boundary condition on the finite points of the boundary. In this
paper, the existence of the approximate solution by the MFS for the Neumann
problems of the modified Helmholtz equation in disk domains is rigorously
demonstrated. We reveal the sufficient condition of the existence of the
approximate solution. Applying Green's theorem to the Neumann problem of the
modified Helmholtz equation, we bound the error between the approximate
solution and exact solution into the difference of the function of the boundary
condition and the normal derivative of the approximate solution by boundary
integrations. Using this estimate of the error, we show the convergence of the
approximate solution by the MFS to the exact solution with exponential order,
that is, $N^2a^N$ order, where $a$ is a positive constant less than one and $N$
is the number of collocation points. Furthermore, it is demonstrated that the
error tends to $0$ in exponential order in the numerical simulations with
increasing number of collocation points $N$.
