Abstract. In the present work, we investigate the applicability of the method of fundamental solutions for the solution of boundary value problems of elliptic partial differential equations and elliptic systems. More specifically, we study whether linear combinations of fundamental solutions can approximate the solutions of the boundary value problems under consideration. In our study, the singularities of the fundamental solutions lie on a prescribed pseudo-boundary -the boundary of a domain which embraces the domain of the problem under consideration. We extend previous density results of Kupradze and Aleksidze, and of Bogomolny, to more general domains and partial differential operators, and with respect to more appropriate norms. Our domains may possess holes and their boundaries are only required to satisfy a rather weak boundary requirement, namely the segment condition. Our density results are with respect to the norms of the spaces C (Ω). Analogous density results are obtainable with respect to Hölder norms. We have studied approximation by fundamental solutions of the Laplacian, biharmonic and m−harmonic and modified Helmholtz and poly-Helmholtz operators. In the case of elliptic systems, we obtain analogous density results for the Cauchy-Navier operator as well as for an operator which arises in the linear theory of thermo-elasticity. We also study alternative formulations of the method of fundamental solutions in cases when linear combinations of fundamental solutions of the equations under consideration are not dense in the solution space. Finally, we show that linear combinations of fundamental solutions of operators of order m ≥ 4, with singularities lying on a prescribed pseudo-boundary, are not in general dense in the corresponding solution space.