I don't think there's anything more conecptually difficult about evaluating <br />
\frac{\langle L(f), f \rangle}{\langle f, f \rangle} where \langle f, g \rangle = \int_a^b f(x)g(x)w(x)\,dx for non-constant w than there is for w = 1.
Given your eigenvalue problem L(f) = \lambda f subject to self-adjoint boundary conditions on f, you approximate f as a linear combination of basis functions which satisfy the boundary conditions: f = \sum_{n = 1}^N a_n \phi_n. Then since L is linear, <br />
L(f) = \sum_{n=1}^N a_n L(\phi_n) = \sum_{n=1}^N a_n\sum_{m=1}^N M_{nm}\phi_m where <br />
L(\phi_n) = \sum_{m=1}^N M_{nm} \phi_m. If you let L(f) = \sum b_m\phi_m(x) then you have now <br />
\mathbf{b}^T = \mathbf{a}^T M and taking the inner product with f you get <br />
\sum_{m} \sum_{n} b_m a_n \langle \phi_m, \phi_n \rangle = \sum_n a_nb_n \|\phi_n\|^2 since the basis functions are orthogonal with respect to \langle \cdot, \cdot \rangle. This gives you the approximation <br />
\lambda = \frac{\langle L(f), f \rangle}{\langle f, f \rangle} \approx<br />
\frac{\mathbf{b}^TD\mathbf{a}}{\mathbf{a}^TD\mathbf{a}}<br />
= \frac{\mathbf{a}^TMD\mathbf{a}}{\mathbf{a}^TD\mathbf{a}} where <br />
D = \operatorname{diag}(\|\phi_1\|^2, \dots, \|\phi_N\|^2). The difficulty with using Bessel functions rather than trigonometric or polynomials is not that the weight function is not constant (there are systems of orthogonal polynomials which also have non-constant weight functions) but in determining <br />
M_{nm} = \langle L(\phi_n), \phi_m \rangle = \int_a^b x L(\phi_n)(x) \phi_m(x)\,dx and \|\phi_n\|^2 = \int_a^b x \phi_n(x)^2\,dx. This is one of the reasons why Bessel functions are less likely to be used in numerical methods than Chebyshev polynomials or finite elements.
If you do want to attempt it, then the relevant chapters of Abramowitz & Stegun are probably a good
place to start.