We present a Matlab implementation of a novel, efficient pseudo-spectral collocation scheme for the solution of stiff, non-local, non-linear PDEs in one and two spatial dimensions. In particular, we compute the non-local terms in real space with the help of a specialised Gauss quadrature. Due to the exponential accuracy of the quadrature and a convenient choice of collocation points near interfaces, we can use grids with a significantly lower number of nodes than most other reported methods. We have demonstrated the capabilities of our numerical methodology by applying it to (Dynamical) Density Functional Theory problems, studying equilibrium and dynamic two-dimensional test cases with single-and multispecies hard-sphere and hard-disc particles modelled with fundamental measure theory, with and without van der Waals attractive forces, in bounded and unbounded physical domains. We have shown that our results satisfy statistical mechanical sum rules. See http://www.research.ed.ac.uk/portal/files/31208790/Pseudospectral_methods_for_density_functional_theory_in_bounded_and_unbounded_domains.pdf for details.See ReadMe.txt