The method of fundamental solutions (MFS) is known as an effective boundary meshless method. However, the formulation of the MFS results in a dense and extremely ill-conditioned matrix. In this paper we investigate the MFS for solving large-scale problems for the nonhomogeneous modified Helmholtz equation. The key idea is to exploit the exponential decay of the fundamental solution of the modified Helmholtz equation, and consider a sparse or diagonal matrix instead of the original dense matrix. Hence, the homogeneous solution can be obtained efficiently and accurately. A standard two-step solution process which consists of evaluating the particular solution and the homogeneous solution is applied. Polyharmonic spline radial basis functions are employed to evaluate the particular solution. Five numerical examples in irregular domains and a large number of boundary collocation points are presented to show the simplicity and effectiveness of our approach for solving large-scale problems.