The equation of radiative transfer is an integro-differential equation in a five-dimensional phase space for the specific intensity of a radiation field. The equation models the transport of the radiation field, the energy loss due to absorption, the energy gain due to emission, and the redistribution of energy due to scattering. In the PN approximation, the specific intensity is replaced by a truncated spherical harmonics expansion, which results in replacing the five-dimensional integro-differential equation by a three-dimensional system of coupled partial differential equations. The resulting system of PDEs is hyperbolic, although the system becomes a parabolic heat equation in the vanishing mean-free path limit (i.e., the scattering dominated regime). A desirable feature of numerical methods for the PN system is that they remain stable and accurate if we fix the mesh parameters and take the vanishing mean-free path limit — in the literature this has been dubbed the “asymptotic-preserving” property. In this work, we develop a Chebyshev pseudo-spectral method for solving the PN system. The time-stepping is done using an L-stable scheme that guarantees that the overall numerical method is asymptotic-preserving. In the multidimensional implementation of the method, we make use of the Radon transform to reduce the computational complexity of the matrix inversion. Several numerical tests are presented in order to demonstrate the feasibility of the resulting method.