Abstract. Many systems of orthogonal polynomials and functions are bases of a variety of function spaces, such as the Hermite and Laguerre functions which are orthogonal bases of L 2 (−∞, ∞) and L 2 (0, ∞), and the Jacobi polynomials which are an orthogonal basis of a weighted L 2 (−1, 1). The associated Legendre functions, and more generally, the spheroidal wave functions are also an orthogonal basis of L 2 (−1, 1).The prolate spheroidal wave functions, which are a special case of the spheroidal wave functions, possess a very surprising and unique property. They are an orthogonal basis of both L 2 (−1, 1) and a subspace of L 2 (−∞, ∞), known as the Paley-Wiener space of bandlimited functions. They also satisfy a discrete orthogonality relation. No other system of classical orthogonal functions is known to possess this strange property. This raises the question of whether there are other systems possessing this property.The aim of the article is to answer this question in the affirmative by providing an algorithm to generate such systems and then demonstrating the algorithm by a new example.