Assume that we have a function $u$ defined in a ball in a discrete way: we know only the values of $u$ in the nodes $(i,j,k)$ of spherical grid, where $i$ is a radius coordinate, $j$ is a coordinate for angle $\varphi$, $k$ is a coordinate for angle $\psi$.
Consider a vector-function $$\nabla u_{i,j,k}=\left(\frac{\partial u}{\partial r}_{i,j,k},\frac{1}{r_i\sin\psi_k}\frac{\partial u}{\partial \varphi}_{i,j,k},\frac{1}{r_i}\frac{\partial u}{\partial \psi}_{i,j,k}\right)-$$ gradient of $u$.
I need to know the values of $\nabla u_{i,j,k}$ on z-axis in cartesian coordinates, which corresponds to $\psi=0$ -- axis in spherical coordinates, but we can not use the formula above, because in case $\psi=0$ the second term turns to infinity.
Actually, we can find the values of $\frac{\partial u}{\partial z}$ with a help of the formula of numerical derivative, but we have a problem with finding $\frac{\partial u}{\partial x}$,$\frac{\partial u}{\partial y}$, because the grid is not rectangular. Could you help me with this stuff and advise me what to do?