Inverse and ill-posed problems which consist of reconstructing the unknown support of a source from a single pair of exterior boundary Cauchy data are investigated. The underlying dependent variable, e.g. potential, temperature or pressure, may satisfy the Laplace, Poisson, Helmholtz or modified Helmholtz partial differential equations (PDEs). For constant coefficients, the solutions of these elliptic PDEs are sought as linear combinations of explicitly available fundamental solutions (free-space Greens functions), as in the method of fundamental solutions (MFS). Prior to this application of the MFS, the free-term inhomogeneity represented by the intensity of the source is removed by the method of particular solutions. The resulting transmission problem then recasts as that of determining the interface between composite materials. In order to ensure a unique solution, the unknown source domain is assumed to be starshaped. This in turn enables its boundary to be parametrised by the radial coordinate, as a function of the polar or, spherical angles. The problem is nonlinear and the numerical solution which minimizes the gap between the measured and the computed data is achieved using the Matlab toolbox routine lsqnonlin which is designed to minimize a sum of squares starting from an initial guess and with no gradient required to be supplied by the user. Simple bounds on the variables can also be prescribed. Since the inverse problem is still ill-posed with respect to small errors in the data and possibly additional ill-conditioning introduced by the spectral feature of the MFS approximation, the least-squares functional which is minimized needs to be augmented with regularizing penalty terms on the MFS coefficients and on the radial function for a stable estimation of these couple of unknowns. Thorough numerical investigations are undertaken for retrieving regular and irregular shapes of the source support from both exact and noisy input data.