I am reading the paper "Discrete Bessel Functions" by R.H. Boyer (Journal of Mathematical Analysis and Applications. Vol 2, Issue 3, June 1961, pg. 509-524) and he begins the paper by discussing Laplace's equation in cylindrical coordinates with axial symmetry: $$\dfrac{\partial}{\partial r} r \dfrac{\partial \phi}{\partial r} + r \dfrac{\partial^2 \phi}{\partial z^2} = 0.$$
He then approximates this PDE with a partial difference equation by using the difference operators $$(\Delta_r \phi)(r,z) = \phi(r+1,z)-\phi(r,z),$$ $$(\nabla_r \phi)(r,z) = \phi(r,z) - \phi(r-1,z),$$ $$(\Delta_z \phi)(r,z) = \phi(r,z+1)-\phi(r,z),$$ and $$(\nabla_z \phi)(r,z) = \phi(r,z) - \phi(r,z-1).$$
For second derivatives, I think he uses the following approximations: $$(\Delta_r \nabla_r \phi)(r,z) = \Delta_r [u(r,z)-u(r-1,z)] = u(r+1,z) - 2u(r,z) + u(r-1,z)$$ and $$(\Delta_z \nabla_z \phi)(r,z) = u(r,z+1) - 2u(r,z) + u(r,z-1)$$
He is considering now to approximate Laplace's equation with these finite difference operators on a grid (I will simplify to a grid with spacing 1 everywhere). His "r" variable becomes "m" and his "z" variable becomes "k". So presumably by plugging these in, he claims that the partial difference equation that results is
$$2 \pi ( m+ \frac{1}{2}) \Delta_r \phi (m,n) - 2 \pi (m - \frac{1}{2}) \Delta_r \phi(m,n) + 2 \pi m(\Delta_z \nabla_z \phi)(m,n)=0.$$
I reckon here that the $\pi$'s out front are related to the axial symmetry, so no problem there. My problem comes from the fact that substituting the difference operators does not lead to this expression. By my own calculation (I can provide details of my derivation if requested), the closest I can get to Boyer's expression after plugging the operators in is $$(m + \frac{1}{2}) \Delta_r \phi(m,n) + \frac{1}{2} \Delta_r \phi(m,n) - m \nabla_r \phi(m,n) + m (\Delta_z \nabla_z \phi)(m,n)=0,$$ which is very close, but there does not appear to me to be a way to resolve it to get his exact expression.
Am I making some kind of error, is there a mistake in Boyer's paper, or is there some convention in dealing with finite differences that I don't know about?
We have $$\left(r \frac{\partial^2}{\partial r^2} +\frac{\partial}{\partial r} +r \frac{\partial^2}{\partial z^2}\right)\phi(r,z) = 0.$$ Letting $h=k=1$, Boyer makes the following natural replacement: $$\begin{eqnarray*} \frac{\partial^2\phi(r,z)}{\partial r^2} &\rightarrow& \phi(m+1,n)-2\phi(m,n)+\phi(m-1,n) \\ \frac{\partial\phi(r,z)}{\partial r} &\rightarrow& \frac{1}{2}(\phi(m+1,n)-\phi(m,n))+\frac{1}{2}(\phi(m,n)-\phi(m-1,n)) \\ &=& \frac{1}{2}(\phi(m+1,n)-\phi(m-1,n)) \\ \frac{\partial^2\phi(r,z)}{\partial z^2} &\rightarrow& \phi(m,n+1)-2\phi(m,n)+\phi(m,n-1). \end{eqnarray*}$$ With this convention we get exactly Boyer's result.