There is one boundary problem
$$\frac{\partial u}{\partial t}= \operatorname{div}\left(a^2 E \nabla u\left(r,\varphi,\psi \right) \right) $$ in a ball $$ B_{1}(0)=\left\{x \in \mathbb{R^3}: \left\| x \right\| \le 1 \right\} $$
(boundary and initial values are now not principal). It's obvious that this equation can be expanded to: \begin{gather} \frac{\partial u}{\partial t}=a^2\left(\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u}{\partial r}\right)+\frac{1}{r\sin\varphi}\frac{\partial }{\partial \varphi}\left(\frac{\sin\varphi}{r} \frac{\partial u}{\partial \varphi} \right)+\frac{1}{r\sin\varphi}\frac{\partial}{\partial \psi}\left(\frac{1}{r\sin\varphi}\frac{\partial u}{\partial \psi} \right)\right).\end{gather} This boundary problem is solved numerically. So, we write an implicit finite-difference scheme
$$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=a^2 \left( \frac{2}{r_{1}+(i-1)\Delta r}\frac{u^{n+1}_{i+1,j,k}-u^{n+1}_{i-1,j,k}}{2\Delta r}+\frac{u^{n+1}_{i+1,j,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i-1,j,k}}{\Delta r^2}\right) +$$ $$+a^2\left(\frac{\cot \left(\varphi_{1}+(j-1)\Delta \varphi \right) ) }{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-u^{n+1}_{i,j-1,k}}{2\Delta \varphi}+\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j-1,k}}{\Delta \varphi^{2}} \right) +$$ $$+a^2\left(\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2 \sin^2\left(\varphi_{1}+(j-1)\Delta \varphi \right)}\frac{u^{n+1}_{i,j,k+1}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j,k-1}}{\Delta \psi^2} \right).$$ The main problem appears in realization of this scheme. The values $r_{1},\varphi_{1}$ -- can be equal to zero in some moments of evaluation ,that's why we face with singularities during evaluation. Could you help me, what should I do with this scheme (perhaps, change it in the moments with singularities) in order to avoid singularities during evaluation?