Let us consider a function $f(r,\theta,\phi)$ that can only be efficiently computed in spherical coordinates. We need its Cartesian gradient. We can implement its analytical spherical derivatives, compute the spherical gradient, $\left(\frac{\partial f}{\partial r},\frac{1}{r}\frac{\partial f}{\partial \theta},\frac{1}{r\sin\theta}\frac{\partial f}{\partial \phi}\right)$, and convert the gradient to Cartesian component.What would you do with the values along the $z$-axis, where $\theta$ is undefined and usually set to 0, and in the origin, where $r$ is 0 as well? More specifically, what would you do with the terms in the denominators $r$ and $r\sin\theta$?
Spherical Gradient along $z$-axis
spherical harmonics