Following is the conversion for spherical to cartesian coordinate \begin{align} x &= r \cos\theta \sin\varphi \\ y &= r \sin\theta \sin\varphi \\ z &= r \cos\varphi \end{align} and we are using the reverse computation to compute spherical coordinate from Cartesian coordinate which is defined as \begin{align} r &= \sqrt{x^2+y^2+z^2} \\ \theta &= \tan^{-1}(y/x) \\ \varphi &= \tan^{-1} \left( \sqrt{x^2+y^2}/z \right) \end{align}
Problem arises when $y$ and $x$ is zero so $\theta$ can take any arbitrary value so during Matlab computations this results in NAN (not a number) which makes $\theta$ discontinuous. Is there any interpolation technique to remove this discontinuity and how to interpret $\theta$ in this case.
$ r^2 = x^2 + y^2 +z^2 $
When $ \varphi = \pi/2 , ( x=y=0, \rightarrow z = \pm r ) $ so the point under reference is either at north pole or at south pole. $\theta$ is indeterminate here. Let not value of $x,y$ go below say $10^{-8} $ or so to determine a $\theta$ value.