I am trying to solve Poisson's equation in an axisymmetric cylindrical domain using finite difference. So I start with my differential equation and boundary conditions and discretize them. However, I'm having trouble thinking of how to discretize the radial portion
\begin{array}{lll} \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial u}{\partial r} \right) + \frac{\partial^2 u}{\partial z^2} = -\frac{\rho}{\varepsilon_0} & \longrightarrow & \frac{r_{i+1/2} u_{i+1, j} - 2 r_{i} u_{i,j} + r_{i-1/2} u_{i-1, j}}{r_i \Delta r^2} + \frac{u_{i, j+1} - 2 u_{i, j} + u_{i, j - 1}}{\Delta z^2} = -\frac{\rho_{i,j}}{\varepsilon_0} \\ \left. \frac{\partial u}{\partial r} \right|_{r = 0} = 0 & \longrightarrow & u_{0,j} = \frac{1}{3} \left(4 u_{1,j} - u_{2,j} \right) ~ \\ \lim_{r\to\infty} u(r, z) = 0 & \longrightarrow & ? ~ \\ u(r, l(r)) = f_3(r) & \longrightarrow & u_{i, l(r_i)/\Delta z} = f_3(r_i) \\ u(r, h) = f_4(r) & \longrightarrow & u_{i, N} = f_4(r_i) \\ \end{array}
where
$$ l(r) = \left\{\begin{aligned} &h - \frac{H}{R} r &&: r \le R\\ &0 &&: r > R \end{aligned} \right.$$ Looking at my old class notes, the professor mentioned a method called irregular singular points as a method to better approximate boundary conditions in (semi-)infinite domains but, I don't understand how I would apply this to 2D systems or to radial systems.
I would appreciate any suggestions.
The actual boundary condition should be
$$\lim_{r \rightarrow \infty} \frac{\partial u}{\partial r} = 0$$
which I approximated at some finite point using the backward second order finite difference.
The actual boundary condition can be seen to be correct by a simple thought experiment by setting $H = 0$ (a rectangular grid) with $\rho = 0$ and we must recover the solution 1D cartesian solution to Laplace's equation.