Abstract-In this paper, we develop a new method for the fast and memory-efficient computation of Slepian functions on the sphere. Slepian functions, which arise as the solution of the Slepian concentration problem on the sphere, have desirable properties for applications where measurements are only available within a spatially limited region on the sphere and/or a function is required to be analyzed over the spatially limited region. Slepian functions are currently not easily computed for large band-limits for an arbitrary spatial region due to high computational and large memory storage requirements. For the special case of a polar cap, the symmetry of the region enables the decomposition of the Slepian concentration problem into smaller subproblems and consequently the efficient computation of Slepian functions for large band-limits. By exploiting the efficient computation of Slepian functions for the polar cap region on the sphere, we develop a formulation, supported by a fast algorithm, for the approximate computation of Slepian functions for an arbitrary spatial region to enable the analysis of modern datasets that support large band-limits. For the proposed algorithm, we carry out accuracy analysis of the approximation, computational complexity analysis, and review of memory storage requirements. We illustrate, through numerical experiments, that the proposed method enables faster computation, and has smaller storage requirements, while allowing for sufficiently accurate computation of the Slepian functions.