frac{\Phi_{i-1,j,k}-2\Phi_{i,j,k}+\Phi_{i+1,j,k}}{(\delta x)^2}, \\
  \Delta_y^2\Phi_{i,j,k} &= \frac{\Phi_{i,j-1,k}-2\Phi_{i,j,k}+\Phi_{i,j+1,k}}{(\delta y)^2}, \\
  \Delta_z^2\Phi_{i,j,k} &= \frac{\Phi_{i,j,k-1}-2\Phi_{i,j,k}+\Phi_{i,j,k+1}}{(\delta z)^2}.
  \label{eq:delz2}\end{aligned}$$

Uniform Cylindrical Grid
------------------------

In uniform cylindrical coordinates, we discretize the computational domain $[R_{\rm min},R_{\rm max}]\times[\phi_{\rm min},\phi_{\rm max}]\times[z_{\rm min},z_{\rm max}]$ with size $L_R\times L_\phi\times L_z$ uniformly into $N_R\times N_\phi\times N_z$ cells. We require that $L_\phi=\phi_{\rm max}-\phi_{\rm min}$ should be an integer fraction of 2$\pi$ to impose periodic boundary condition along the azimuthal direction. We define the face-centered radial and azimuthal coordinates as $R_{i+1/2} = R_{\rm min} + i \delta R$ and $\phi_{j+1/2} = \phi_{\rm min} + j \delta \phi$, where $\delta R \equiv L_R/N_R$ and $\delta \phi \equiv L_\phi/N_\phi$. Unlike in Cartesian coordinates, the definition of the radial cell-center is ambiguous in cylindrical coordinates because the geometric center does not coincide with the volumetric center. When the radial grid is uniform, finite difference of quantities defined at the geometric centers can retain second-order accuracy. We thus define the cell-centered coordinates as $R_{i} = (R_{i-1/2} + R_{i+1/2})/2$ with $i=1,2,\cdots,N_R$ and $\phi_{j} = (\phi_{j-1/2} + \phi_{j+1/2})/2$ with $j=1,2,\cdots,N_\phi$. Discretization in the vertical direction is the same as in the uniform Cartesian coordinates.

The second-order finite-difference approximation to Equation can be written as $$\label{eq:fd_uniform_cylindrical}
  \left( \Delta_R^2 + \Delta_\phi^2 + \Delta_z^2 \right)\Phi_{i,j,k} = 4\pi G \rho_{i,j,k},$$ where the difference operators $\Delta_R^2$ and $\Delta_\phi^2$ are defined by $$\begin{aligned}
  \Delta_R^2\Phi_{i,j,k} =& \frac{\Phi_{i-1,j,k}-2\Phi_{i,j,k}+\Phi_{i+1,j,k}}{(\delta R)^2}
    + \frac{\Phi_{i+1,j,k}-\Phi_{i-1,j,k}}{2R_i\de