I'm trying to solve a heat equation in cylindrical coordinates
$$\dfrac{\partial u}{\partial t} = a \left(\dfrac{\partial^2 u}{\partial r^2} + \dfrac{1}{r} \dfrac{\partial u}{\partial r} + \dfrac{1}{r^2} \dfrac{\partial^2 u}{\partial \phi^2} + \dfrac{\partial^2 u}{\partial z^2} \right),~~~~~~~~~~~~~~ (r, \phi, z) \in \Omega.$$ with Neumann boundary conditions $$ \left.\dfrac{\partial u}{\partial n} \right|_{\partial \Omega} = f(u, r, \phi, z),$$ and initial conditions $$u(0, r, \phi, z) = g(r, \phi, z).$$ The region $\Omega$ is just a cylinder $$\Omega =\{ (r, \phi, z) : 0\leq r \leq R, ~~ 0 \leq \phi < 2\pi, ~~ 0 \leq z \leq Z\}.$$ Functions $f$ and $g$ are known and $f$ is a non-linear function with respect to $u$. Due to its non-linearity I cannot solve the problem analytically and trying to use a finite-difference method: $$\begin{array}{lcl} t_\tau = \Delta t \cdot \tau, & \tau=0,\cdots, T, &\\ r_i = \Delta r \cdot i, & i=0,\cdots, I, & (r_I = R),\\ \phi_j = \Delta \phi \cdot j, & i=0,\cdots, J, & (\phi_J = 2\pi - \Delta \phi),\\ z_k = \Delta z \cdot k, & k=0,\cdots, K, & (z_K = Z), \end{array}$$ And $$u^\tau_{i,j,k} = u(t_\tau, r_i, \phi_j, z_k).$$ So $$\dfrac{u^{\tau+1}_{i,j,k} - u^{\tau}_{i,j,k}}{\Delta t} = a\left( \dfrac{u^\tau_{i+1,j,k} - 2 u^\tau_{i,j,k} + u^\tau_{i-1,j,k}}{\Delta r^2} + \dfrac{1}{\Delta r \cdot i} \dfrac{u^\tau_{i+1,j,k} - u^\tau_{i-1,j,k}}{2 \Delta r} + \dfrac{1}{(\Delta r \cdot i)^2} \dfrac{u^\tau_{i,j+1,k} - 2 u^\tau_{i,j,k} + u^\tau_{i,j-1,k}}{\Delta \phi^2} + \dfrac{u^\tau_{i,j,k+1} - 2 u^\tau_{i,j,k} + u^\tau_{i,j,k-1}}{\Delta z^2}\right).$$
Next, I use boundary conditions and periodicity of the solution to find $u^\tau_{I+1, j, k}$, $~~u^\tau_{-1, j, k}, ~~$ etc. For example: $$ \dfrac{u^\tau_{I+1,j,k} - u^\tau_{I-1,j,k}}{2 \Delta r} = f(u^\tau_{I,j,k}, \Delta r \cdot I, \Delta \phi \cdot j, \Delta z \cdot k).$$
However, I have a hard time with Laplacian at $r=0$, i.e., writing the difference equations for $u^\tau_{0, j, k}$.
I know that if the problem was symmetric with respect to $\phi$, then $$\dfrac{\partial^2 u}{\partial \phi^2} \equiv 0.$$ Moreover, $$\left. \dfrac{\partial u}{\partial r}\right|_{r=0} = 0,$$ and, therefore, $$\left. \dfrac{1}{r} \dfrac{\partial u}{\partial r}\right|_{r=0} = \dfrac{\partial^2 u}{\partial r^2}.$$ So, the heat equation at $r=0$ would transform into $$\dfrac{\partial u}{\partial t} = a \left(2\dfrac{\partial^2 u}{\partial r^2} + \dfrac{\partial^2 u}{\partial z^2} \right),$$ which would allow me to just use a separate equations for $u^\tau_{0, j, k}$.
What should be the heat equation at $r=0$ in my more general case?
I think I came up with the answer myself. I'll leave it here if anyone meets the similar problem.
Suppose $f(x, y, z)$ is a sufficiently smooth function defined in Cartesian coordinates and $g(r, \phi, z) = f(r \cos\phi, r\sin\phi, z)$ is the same function in cylindrical coordinates. We have $$\left. \nabla^2 f \right|_{\begin{array}{l}x=r\cos\phi\\y=r\sin\phi\\z=z\end{array}} = \nabla^2 g.$$ Using second-order central difference approximation for the second derivatives: $$\left. \nabla^2 f \right|_{\begin{array}{l}x=0\\y=0\end{array}} \approx \dfrac{f(\Delta x, 0, z) - 2f(0, 0, z) + f(-\Delta x, 0, z)}{(\Delta x)^2} + \dfrac{f(0, \Delta y, z) - 2f(0, 0, z) + f(0, -\Delta y, z)}{(\Delta y)^2} + \dfrac{f(0, 0, z+\Delta z) - 2f(0, 0, z) + f(0, 0, z-\Delta z)}{(\Delta z)^2} = \dfrac{g(\Delta r, 0, z) - 2g(0,\phi,z)+g(\Delta r,\pi,z)}{(\Delta r)^2} + \dfrac{g(\Delta r, \frac{\pi}{2}, z) - 2g(0,\phi,z)+g(\Delta r,\frac{3\pi}{2},z)}{(\Delta r)^2} + \dfrac{g(0, 0, z+\Delta z) - 2g(0, \phi, z) + g(0, 0, z-\Delta z)}{(\Delta z)^2},$$ where $\phi$ can be anything, since $g(0,\phi_1,z)=g(0,\phi_2,z)$ for any $\phi_1, \phi_2$ as they correspond to the same point.
Thus, the finite difference equation that corresponds to the heat equation in the cylindrical coordinates at $r=0$ is the following (using notations from the question): $$ \dfrac{u^{\tau+1}_{0,j,k} - u^{\tau}_{0,j,k}}{\Delta t} = a \left( \dfrac{u^\tau_{1,0,k} - 2u^\tau_{0,j,k} + u^\tau_{1,\beta,k}}{(\Delta r)^2} + \dfrac{u^\tau_{1,\alpha,k} - 2u^\tau_{0,j,k} + u^\tau_{1,\gamma,k}}{(\Delta r)^2} + \dfrac{u^\tau_{0,j,k+1} - 2u^\tau_{0,j,k} + u^\tau_{0,j,k-1}}{(\Delta z)^2}\right),$$ where $\alpha \cdot \Delta\phi = \pi/2$, $~~\beta \cdot \Delta\phi = \pi,~~$ and $~~\gamma \cdot \Delta\phi = 3\pi/2~~$ (so the grid should be appropriately chosen). And once again, $u^\tau_{0,j,k}$ doesn't really depend on $j$.