An efficient and accurate method, based on the weighted averages (WA) extrapolation technique, is presented for the evaluation of semi-infinite range integrals involving products of Bessel functions of arbitrary order. The method requires splitting the integration interval into a finite and an infinite part. The integral over the first finite part is computed using an adaptive quadrature rule based on Patterson formulas. For the evaluation of the remaining integral, the strongly irregular oscillatory behavior of the product of two Bessel functions is first represented as a sum of two asymptotically simply oscillating functions. Then, by applying the integration-then-summation technique, a sequence of partial integrals is obtained, and its convergence is accelerated with the help of WA. Details and possible complications involved in the method are addressed. Finally, the excellent performance of the proposed method is demonstrated throughout several numerical examples.