In this paper, we study the Helmholtz equation by the method of fundamental solutions (MFS) using Bessel and Neumann functions. The bounds of errors are derived for bounded simply-connected domains, while the bounds of condition number are derived only for disk domains. The MFS using Bessel functions is more efficient than the MFS using Neumann functions. Note that by using Bessel functions, the radius R of the source nodes is not necessarily to be larger than the maximal radius rmax of the solution domain. This is against the well-known rule: rmax < R for the MFS. Numerical experiments are carried out, to support the analysis and conclusions made. This is the first novelty in this paper. The error analysis for the Helmholtz equation is more complicated than that for the modified Helmholtz equation in [Appl. Math. Comp., 305(2017):262-281], since the Bessel functions Jn(x) have infinite zeros. We consider the curial and degenerate cases when Jn(kR) ≈ 0 and Jn(kρ) ≈ 0. There exist few reports for the analysis for such a degeneracy. The error bounds are also explored for bounded simply-connected domains. The second novelty of this paper is for the analysis of the MFS in degeneracy. For the MFS using Neumann functions, the rule of the MFS, rmax < R, must obey. This paper is the first time to discover that the MFS using Bessel and Neumann functions suffer from the spurious eigenvalues. The spurious eigenvalues are not the true eigenvalues of the corresponding eigenvalue problems, but the correct solutions can not be obtained due to either algorithm singularity or divergence of numerical solutions. For the method of particular solutions (MPS), however, the source nodes disappear. In this paper, we will briefly provide the analysis of the MFS using Neumann functions, and the polynomial convergence can be achieved for bounded simple-connected domains. The analysis of the MFS using Neumann functions and numerical comparisons for different methods are the third contribution in this paper.