We construct solutions of the paraxial and Helmholtz equations that are polynomials in their spatial variables. These are derived explicitly by using the angular spectrum method and generating functions. Paraxial polynomials have the form of homogeneous Hermite and Laguerre polynomials in Cartesian and cylindrical coordinates, respectively, analogous to heat polynomials for the diffusion equation. Nonparaxial polynomials are found by substituting monomials in the propagation variable z with reverse Bessel polynomials. These explicit analytic forms give insight into the mathematical structure of paraxially and nonparaxially propagating beams, especially in regard to the divergence of nonparaxial analogs to familiar paraxial beams.