I have the following Sturm-Liouville eigenvalue problem,
$$\frac{d}{dr}\left(r \frac{d u}{dr} \right) + k^2 r u = 0, \, \, u(R) = 0, \, u(0) < \infty$$
where $k^2$ the eigenvalues and $u(r)$ the associated eigenfunctions. I was able to prove that the operator $L = \frac{d}{dr}\left(r \frac{d }{dr} \right), (Lu + k^2ru=0),$ is autoadjoint under the scalar product $\int_0^R f\, g dr$, and that the eigenfunctions of different eigenvalues are orthogonal under $\int_0^r f \, g \, r dr$.
However I want to solve the eigenvalue problem, so I want to find the eigenvalues and its corresponding eigenfunctions. This is actually the Bessel equation of zero order, so I have tried to solve this with a Frobenius series; $u(r) = \sum_{l = 0}^{\infty} a_l z^{l+c}$ (due to we have a singular regular point at $r=0$), so finally I have obtained the following expression,
$u(r) = \sum_{l = 0}^{\infty} \frac{(-1)^l}{(l!)^2} \left( \frac{k r}{2}\right)^{2l}$
How do I obtain the eigenvalues $k^2$?
The eigenvalues are determined by the boundary condition $u(R) = 0$. As you have obtained $u(r) = J_0(k r)$, these are the values for which the Bessel function $J_0(k R)$ is zero. The zeroes of the Bessel function cannot be expressed in elementary functions. However, you can use the asymptotic approximations of $J_0(z)$ for small and large values of $z$ to get asymptotic expressions of the eigenvalues (see eg. here). For more information on zeroes of $J_0(z)$, see DLMF.