Inhomogeneous elliptical inclusions with partial differential equations have aroused appreciable concern in many disciplines. In this paper, a pseudo-spectral collocation method, based on Fourier basis functions, is proposed for the numerical solutions of two- (2D) and three-dimensional (3D) inhomogeneous elliptic boundary value problems. We describe how one can improve the numerical accuracy by making some extra “reconstruction techniques” before applying the traditional Fourier series approximation. After the particular solutions have been obtained, the resulting homogeneous equation can then be calculated using various boundary-type methods, such as the method of fundamental solutions (MFS). Using Fourier basis functions, one does not need to use large matrices, making accrual computations relatively fast. Three benchmark numerical examples involving Poisson, Helmholtz, and modified-Helmholtz equations are presented to illustrate the applicability and accuracy of the proposed method.