A single closed expression is derived for the product-sum conversion of spherical harmonics and their derivatives to any order. The equations for thermal convection in a sphere at an infinite Prandtl number are transformed to a set of quadratic first-order ordinary differential equations suitable for numerical integration with respect to radius. The product-sum formulas are applied to the quadratic terms arising from advection, dissipation, and the linear temperature dependence of viscoity.