We present an algorithm to compute integrals of the form ∫ ∞ 0 x m ∏ k i = 1 J ν i ( a i x ) dx with J ν i ( x ) the Bessel function of the first kind and (real) order ν i . The parameter m is a real number such that ∑ i ν i + m > −1 and the coefficients a i are strictly positive real numbers. The main ingredients in this algorithm are the well-known asymptotic expansion for J ν i ( x ) and the observation that the infinite part of the integral can be approximated using the incomplete Gamma function Γ( a , z ). Accurate error estimates are included in the algorithm, which is implemented as a MATLAB program.