To make it easier to understand, you’ll want to write the first two terms together in a more symmetrical form:
\frac{\partial^2 u}{\partial r^2} + \frac{1}{r} \frac{\partial u}{\partial r} = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r} \right)
which is now clearly a symmetric operator under the weighted inner product \langle u, v \rangle = \int u v \, r dr with appropriate boundary conditions. (Understanding the symmetry of the infinite-dimensional operator is really critical to discretizing it properly.)
Now, suppose you discretize the \partial / \partial r first-derivative operator by an (n+1)\times n matrix D which takes in n values \vec{u} on a grid G, plus two points on the boundaries (I’m assuming Dirichlet boundaries for simplicity) and returns the centered differences at n+1 points G' (halfway between the original grid points including the boundaries). After a little algebra, the discretized second-derivative operator above is then the n\times n matrix operation:
\left. \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r} \right) \right|_G \approx -R^{-1} D^T R' D \vec{u}
where R is a diagonal n \times n matrix of r values on G, R' is a diagonal (n+1) \times (n+1) matrix of r values on the grid G', and -D^T turns out to be the center-difference operator on G' that gives you the approximate derivative on G. But -R^{-1} D^T R' D is clearly similar to a symmetric (and positive-definite!) matrix under a diagonal change of basis \vec{v} = R^{1/2}\vec{u}, which corresponds exaclty to the r Jacobian factor in the inner product for the continuous operator. Equivalently, it’s symmetric under a weighted inner product \langle \vec{u}, \vec{v} \rangle = \vec{u}^T R \vec{v}. (Any Krylov method, such as conjugate gradient, is easily adapted to support modified inner products.)
The above works if your domain is an annulus with Dirichlet boundaries, but if your domain includes the origin, you need to be a little more careful in setting up finite differences because of the coordinate singularity. You typically want to impose Neumann and not Dirichlet boundaries at the origin, and you may want to choose your discretization so that the grid G does not include the origin (so that you don’t divide by zero), or otherwise you can work out the derivative rule at the origin by a limiting procedure.
(This can be combined with a periodic center-difference discretization in θ, which is trivial symmetrical, using a Kronecker product, and the Kronecker product of symmetric matrices is symmetric.)
