The long-standing problem of implementing the PN method effectively for spherical geometry is revisited in this work. It is shown that a least-squares approach to the method resolves to a great extent the numerical instability reported for the first time by Aronson in 1984. In the proposed version of the method, a small loss of accuracy is still observed for intermediate orders of the approximation, but in high order (typically N ≥ 199), full accuracy is recovered, and the method can be used with confidence even for extremely high orders of the approximation. Numerical results of benchmark quality are tabulated for the quantities of interest for two basic transport problems in spherical geometry: the albedo problem for a sphere and the critical-sphere problem, both including cases that show the effects of scattering anisotropy described by the binomial law.