I recommend you to visit https://dlmf.nist.gov which includes a lot of information about different OP. In particular Bessel functions... https://dlmf.nist.gov/10
When r_1 and r_2 are different, the integral is calculated explicitly in terms of a sum of products of simple rational functions and values of J_{n+1/2} and J_{n-1/2} at the points ar_1, ar_2,br1 and br2, where J_n denote the Bessel function of the first kind. The idea is straightforward. First of all, having in mind the simple relation between j_n and J_n, you need to integrate the product J_n(r_1 x) J_n(r_2 x) against the weight x and not against x^2 as in your original problem. Now look at http://physics.ucsc.edu/~peter/116C/bess_orthog.pdf to understand how the orthogonality of the Bessel functions is obtained. Summarising, when you obtain (10), with a and b there being your r_1 and r_2, you have to integrate that relation in [a,b]. Finally, you have to use some simple formulae for the derivative of the Bessel function to obtain the following result, which is given also via mathematica:
-(a^2 π (r2 BesselJ[-(1/2) + n, a r2] BesselJ[1/2 + n, a r1] - r1 BesselJ[-(1/2) + n, a r1] BesselJ[1/2 + n, a r2]))/(2 Sqrt[a r1] Sqrt[ a r2] (r1^2 - r2^2)) + (b^2 π (r2 BesselJ[-(1/2) + n, b r2] BesselJ[1/2 + n, b r1] - r1 BesselJ[-(1/2) + n, b r1] BesselJ[1/2 + n, b r2]))/(2 Sqrt[b r1] Sqrt[b r2] (r1^2 - r2^2))
Obviously, when [a,b]=[0,1], the latter is equal to zero, which is exactly the orthogonality of the Bessel functions.