A new set of polynomial functions that can be used in spectral expansions of C∞ functions in polar coordinates (r, φ) is defined by a singular Sturm-Liouville equation. With the use of the basis functions, the spectral representations remain analytic at the pole despite the coordinate singularity because the pole condition is exactly satisfied at the origin for ag azimuthal modes, not just a few of the gravest modes (which is the usual case). Based on recurrence relations, fast and stable numerical operators for 1/r, r(d/dr), the Laplacian and Helmholtz operators and their inverses are developed. Although the spacings in the azimuthal direction of the collocation points near the origin are small (i.e., α 1/M2, where M is the number of radial modes), the explicit numerical method for Euler's equation is not stiff at the origin. Namely, the CFL number σ is O(1) where the grid size in is defined as π/M (i.e., the maximum allowable timestep is proportional to 1/M, not 1/M2).