Proving that $\int_0^\infty\frac{J_{2a}(2x)~J_{2b}(2x)}{x^{2n+1}}~dx~=~\frac12\cdot\frac{(a+b-n-1)!~(2n)!}{(n+a+b)!~(n+a-b)!~(n-a+b)!}$
Using the product cited above, viz.
\begin{align} J_{\mu}(x) \, J_{\nu}(x) &= \sum_{n=0}^{\infty} \frac{(-1)^{n} \, \Gamma(2n+\mu+\nu+1) \, \left(\frac{x}{2}\right)^{2n+\mu+\nu}}{n! \, \Gamma(\mu+\nu+n+1) \, \Gamma(\mu+n+1) \, \Gamma(\nu+n+1)} \\ \implies \frac{J_{\mu}(\sqrt x) \, J_{\nu}(\sqrt x)}{x^{\mu/2+\nu/2}} &= \sum_{n=0}^{\infty} \frac{(-1)^{n} \, \Gamma(2n+\mu+\nu+1) \, x^n}{2^{2n+\mu+\nu} n! \, \Gamma(\mu+\nu+n+1) \, \Gamma(\mu+n+1) \, \Gamma(\nu+n+1)} \\ \end{align}
we can invoke Ramanujan's Master Theorem, to get
\begin{align} \int_0^\infty \frac{J_{\mu}(\sqrt x) \, J_{\nu}(\sqrt x)}{x^{\mu/2+\nu/2}} x^{s-1}\,dx = \Gamma(s) \frac{\Gamma(-2s+\mu+\nu+1)}{2^{-2s+\mu+\nu}\, \Gamma(\mu+\nu-s+1) \, \Gamma(\mu-s+1) \, \Gamma(\nu-s+1)} \end{align}
If $x=4 u^2 \implies dx = 8u\,du$, the LHS becomes
\begin{align} 2^{2s+1-\mu-\nu}\int_0^\infty J_{\mu}(2u) \, J_{\nu}(2u) u^{2s-\mu-\nu-1}\,du \end{align}
Therefore,
\begin{align} \int_0^\infty J_{\mu}(2u) \, J_{\nu}(2u) u^{2s-\mu-\nu-1}\,du = \frac{\Gamma(s) \Gamma(-2s+\mu+\nu+1)}{2\, \Gamma(\mu+\nu-s+1) \, \Gamma(\mu-s+1) \, \Gamma(\nu-s+1)} \\ \int_0^\infty J_{2\mu}(2u) \, J_{2\nu}(2u) u^{2s-2\mu-2\nu-1}\,du = \frac{\Gamma(s) \Gamma(-2s+2\mu+2\nu+1)}{2\, \Gamma(2\mu+2\nu-s+1) \, \Gamma(2\mu-s+1) \, \Gamma(2\nu-s+1)} \end{align}
Set $s=\mu+\nu-n$.
\begin{align} \int_0^\infty \frac{J_{2\mu}(2u) \, J_{2\nu}(2u)}{u^{2n+1}}\,du &= \frac{\Gamma(\mu+\nu-n) \Gamma(2n+1)}{2\, \Gamma(\mu+\nu+n+1) \, \Gamma(-\mu+\nu+n+1) \, \Gamma(\mu-\nu+n+1)} \\ &= \frac{(\mu+\nu-n-1)!(2n)!}{2(\mu+\nu+n)!(-\mu+\nu+n)!(\mu-\nu+n)!} \end{align}
This is your solution.
To save from typing a result a highlight will be listed for now.
- The integral in question is a reduction of the more general integral \begin{align} \int_{0}^{\infty} \frac{J_{\mu}(at) \, J_{\nu}(bt)}{t^{\lambda}} \, dt = \frac{b^{\nu} \Gamma\left( \frac{\mu + \nu - \lambda +1}{2}\right)}{2^{\lambda} \, a^{\nu - \lambda +1} \, \Gamma\left( \frac{\lambda + \mu - \nu +1}{2} \right) } \, {}_{2}F_{1}\left( \frac{\mu+\nu-\lambda+1}{2}, \frac{\nu-\lambda-\mu+1}{2}; \nu+1; \frac{b^{2}}{a^{2}} \right). \end{align} When $\mu \rightarrow 2 \mu$, $\nu \rightarrow 2 \nu$, $a=b=2$, $\lambda = 2n+1$ the result is obtained
- A method to obtain the above listed result is to consider the integral as \begin{align} \int_{0}^{\infty} \frac{J_{\mu}(at) \, J_{\nu}(bt)}{t^{\lambda}} \, dt = \lim_{s \rightarrow 0} \, \int_{0}^{\infty} e^{-s t} \, t^{-\lambda} \, J_{\mu}(at) \, J_{\nu}(bt) \, dt \end{align}
- See G. N. Watson's Bessel function Book section 13.4, p.401
Edit: The product of two Bessel functions, as required by this problem, is \begin{align} J_{\mu}(x) \, J_{\nu}(x) = \sum_{n=0}^{\infty} \frac{(-1)^{n} \, \Gamma(2n+\mu+\nu+1) \, \left(\frac{x}{2}\right)^{2n+\mu+\nu}}{n! \, \Gamma(\mu+\nu+n+1) \, \Gamma(\mu+n+1) \, \Gamma(\nu+n+1)} \end{align} When $x = 2t$ it is seen that \begin{align} \int_{0}^{\infty} e^{-st} \, t^{- \lambda} \, J_{\mu}(2t) \, J_{\nu}(2t) \, dx = \sum_{n=0}^{\infty} \frac{(-1)^{n} \, \Gamma(2n+\mu+\nu+1) \, \Gamma(2n+\mu+\nu-\lambda+1) \, \left(\frac{1}{s}\right)^{2n+\mu+\nu-\lambda+1}}{n! \, \Gamma(\mu+\nu+n+1) \, \Gamma(\mu+n+1) \, \Gamma(\nu+n+1)} \end{align} Reducing this series and a possible transformation of the resulting hypergeometric series along with the limiting value for $s$ will yield the desired result.
Note: Here is an outline of an alternate proof, which could be helpful.
In section $21$ in Bessel Functions and Their Applications B.G. Korenev introduces the so-called discontinuous integrals of Weber-Sonin-Schafheitlin \begin{align*} \int_0^{\infty}\frac{J_\mu(at)J_\nu(bt)}{t^\lambda}\,dt \end{align*} He notes, that some special cases of these integrals were calculated by Weber; for all values of the parameters $\lambda,\mu,\nu$ such that this integral converges it has been calculated by Sonin; Schafheitlin in his investigations has noted that at $a=b$ this integral has a discontinuity.
Korenev provides in his book a summary of the main points of the proof and refers at the end of the section to volume II of H. Bateman's classic Higher Transcendental Function for more details.
The relevant part in volume II is section $7.7.4$ The discontinuous integral of Weber and Schafheitlin. We find in
[H. Bateman] Volume II, Section $7.7.4$:
We shall now investigate the integral \begin{align*} \int_0^{\infty}\frac{J_\mu(at)J_\nu(bt)}{t^\rho}\,dt \end{align*} in which are $a,b$ are positive real. It turns out that even when the integral converges for all positive $a$ and $b$, its analytic expression is different, according as $a$ is smaller, equal to, or larger than $b$. The results are
\begin{align*} \int_0^\infty&\frac{J_\mu(at)J_\nu(bt)}{t^\rho}\,dt\\ &=\frac{a^{\mu}\Gamma\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\mu-\frac{1}{2}\rho\right)} {2^\rho b^{\mu-\rho+1}\Gamma(\mu+1) \Gamma\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\rho-\frac{1}{2}\mu\right)}\tag{1}\\ &\quad\cdot {}_{2}F{_1}\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\mu-\frac{1}{2}\rho,\frac{1}{2}+\frac{1}{2}\mu-\frac{1}{2}\nu-\frac{1}{2}\rho;\mu+1;\frac{a^2}{b^2}\right)\\ &\\ &\Re(\nu+\mu-\lambda+1)>0,\quad \Re(\rho)>-1,\quad0<a<b \end{align*}
with a corresponding expression for $0<b<a$ [interchange $a$ and $b$] and
\begin{align*} \int_0^\infty&\frac{J_\mu(at)J_\nu(at)}{t^\rho}\,dt\tag{2}\\ &=\frac{{\left(\frac{1}{2}a\right)}^{\rho-1} \Gamma(\rho) \Gamma\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\mu-\frac{1}{2}\rho\right) }{2 \Gamma\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\rho-\frac{1}{2}\mu\right) \Gamma\left(\frac{1}{2}+\frac{1}{2}\nu+\frac{1}{2}\rho+\frac{1}{2}\mu\right) \Gamma\left(\frac{1}{2}-\frac{1}{2}\nu+\frac{1}{2}\rho+\frac{1}{2}\mu\right) }\\ &\\ &\qquad\qquad\qquad\Re(\nu+\mu+1)>\Re(\rho)>0,\quad a>0 \end{align*}
Note, that substituting in (2) \begin{align*} a\rightarrow 2,\quad t\rightarrow x,\quad \mu\rightarrow2a,\quad \nu\rightarrow2b,\quad \rho\rightarrow2n+1 \end{align*} we obtain \begin{align*} \int_0^\infty&\frac{J_{2a}(2x)J_{2b}(2x)}{x^{2n+1}}\,dx\tag{3}\\ &=\frac{1}{2}\frac{ \Gamma(2n+1) \Gamma\left(b+a-n\right) }{ \Gamma\left(b+n-a+1\right) \Gamma\left(b+n+a+1\right) \Gamma\left(-b+n+a+1\right) } \end{align*} corresponding to OPs expression.
In order to proof the formulae above Bates continues as following
([Bates] adapted:) We use in the integrand of (1) a generalization of Neumann's formula, \begin{align*} J_{\mu}(at)J_{\nu}(bt)=\frac{1}{\pi}(2t)^{\mu+\nu}&a^\mu b^\nu \int_{-\frac{1}{2}\pi}^{\frac{1}{2}\pi}e^{i\theta(\mu-\nu)}\left(\frac{\cos \theta}{\lambda t}\right)^{\mu+\nu} J_{\mu+\nu}(\lambda t) \,d\theta\tag{4}\\ &\\ \Re(\mu+\nu)>-1,&\quad\lambda=\left(2\cos\theta\left(a^2e^{i\theta}+b^2e^{-i\theta}\right)\right)^\frac{1}{2} \end{align*} interchange the order of integration, evaluate the integral with respect to $t$ by \begin{align*} \int_{0}^{\infty}J_\mu(at)t^{\rho -1}\,dt&= 2^{\rho -1}a^{-\rho}\frac{\Gamma\left(\frac{1}{2}\mu+\frac{1}{2}\rho\right)}{\Gamma\left(1+\frac{1}{2}\mu-\frac{1}{2}\rho\right)} \tag{5}\\ &\\ -\Re (\mu)&<\Re (\rho)<3/2,\quad a>0 \end{align*} and obtain \begin{align*} \int_0^{\infty}&\frac{J_{\mu}(at)J_{\nu}(bt)}{t^{\rho}}\,dt= \frac{1}{\pi}a^{\mu}b^{\nu}2^{\frac{1}{2}\mu+\frac{1}{2}\nu-\frac{1}{2}\rho-\frac{1}{2}} \frac{\Gamma\left(\frac{1}{2}\mu+\frac{1}{2}\nu-\frac{1}{2}\rho+\frac{1}{2}\right)} {\Gamma\left(\frac{1}{2}\mu+\frac{1}{2}\nu+\frac{1}{2}\rho+\frac{1}{2}\right)}\\ &\qquad\cdot\int_{-\frac{1}{2}\pi}^{\frac{1}{2}\pi}e^{i\theta(\mu-\nu)}\left(\cos \theta\right)^{\frac{1}{2}(\mu+\nu+\rho-1)} \left(a^2e^{i\theta}+b^2e^{-i\theta}\right)^{\frac{1}{2}(\rho-\mu-\nu-1)}\,d\theta \end{align*} But the integrand on the right-hand side is expressible as a hypergeometric function ${}_2F_1$, compare volume I, 2.4(11) and we immediately obtain the expressions (1) and (2) according as $b>a$ or $b=a$.
Please note, it is more an outline of a proof since some details are omitted. But it's hopefully useful for further study.
Hint: When looking for an answer, I also found in chapter $4$, Bessel Functions and Confluent Hypergeometric Functions in Special Functions by J. Andrews, etal.
Exercise 15: Show that when $a=b$ in Exercise $14$, the result is \begin{align*} \int_0^{\infty}\frac{J_{\alpha-\beta}(at)-J_{\gamma-1}(at)}{t^{\gamma-\alpha-\beta}}\,dt &=\frac{(a/2)^{\gamma-\alpha-\beta-1}\Gamma(\gamma-\alpha-\beta)\Gamma(\alpha)} {2\Gamma(1-\beta)\Gamma(\gamma-\alpha)\Gamma(\gamma-\beta)} \end{align*} provided that $\Re (\alpha)>0$ and $\Re( \gamma-\alpha-\beta) > 0$.
In Exercise $14$ which is a more general variant of the integral above there is a also reference pointing to Watson [1944] section 13.4, as we already now from the answer of @Leucippus. In fact the expression in Exercise $15$ is stated as formula (1) on p. $403$ in Watsons book.
So, it could be worthwile to read chapter $4$ about Bessel functions in J. Andrews's book. It consists of roughly $40$ pages, which provides all knowledge necessary to fully grasp Watson's proof.