Pseudospectral methods utilizing prolate spheroidal wave functions as basis functions have been shown to possess advantages over the conventional pseudospectral methods based on trigonometric and orthogonal polynomials. However, the spectral differentiation and interpolation steps of the prolate pseudospectral method involve matrix-vector products, which, if evaluated directly, entail O(N2) memory requirement and computational complexity (where N is the number of unknowns utilized for discretization and interpolation). In this work we show that the fast multipole method (FMM) can be used to reduce the memory requirement and computational complexity of the prolate pseudospectral method to O(N). Example simulation results demonstrate the speed and accuracy of the resulting fast prolate pseudospectral solver. © 2006 Society for Industrial and Applied Mathematics.