An integral that comes up when discussing the orthogonality relations of the Bessel function of the first kind is $$ \int x J_p(\alpha x ) J_p(\beta x) \, \mathrm dx = \frac{x\left(\alpha J'_p(ax)J_p(\beta x) - \beta J_p(ax) J'_p(\beta x) \right)}{\beta^2 - \alpha^2} + C $$
The proof, which uses the defining differential equation of Bessel function of the first kind, is quite simple.
I would like to confirm that the derivative of the right side of the equation is indeed $$x J_p(\alpha x ) J_p(\beta x)$$
I tried using identities like $$J_{p-1}(x) - J_{p+1}(x) = 2 J'_p(x)$$ and $$J_{p-1}(x) + J_{p+1}(x) = \frac{2p}{x} J_p(x), $$ but I feel that I'm just going in circles.
This s not actually anything to do with the relations between Bessel functions with different $p$: it arises purely from the fact that Bessel functions solve Bessel's equation. This can be written in Sturm–Liouville form as $$ -(xu')' + \frac{p^2}{x}u=x u. $$ Differentiating the right-hand side of the given expression using the product rule judiciously gives \begin{align} &(\alpha x J_p'(\alpha x)J_p(\beta x)-\beta x J_p'(\beta x)J_p(\alpha x))' \\ &= \alpha( x J_p'(\alpha x))'J_p(\beta x) + \alpha x J_p'(\alpha x) (J_p(\beta x))' - \beta (x J_p'(\beta x))' J_p(\alpha x)-\beta x J_p'(\beta x)(J_p(\alpha x))' \\ &= \alpha (x J_p'(\alpha x))'J_p(\beta x) - \beta (x J_p'(\beta x))' J_p(\alpha x) + \alpha\beta x J_p'(\alpha x) J_p'(\beta x) -\alpha\beta x J_p'(\beta x)J_p'(\alpha x) \\ &= ( x J_p(\alpha x)')'J_p(\beta x) - (x J_p(\beta x)')' J_p(\alpha x) \end{align} Now, if we change variables in the Sturm–Liouville equation to $x=ay$, then $ u' = \frac{1}{a} \frac{du}{dy} $, so $u(ay)$ satisfies $$ -\frac{1}{a}\frac{d}{dy}\left(y\frac{du}{dy}\right) + \frac{p^2}{ay}u =ay u, $$ or $$ -\frac{d}{dy}\left(y\frac{du}{dy}\right) + \frac{p^2}{y}u=a^2y u, \tag{*} $$ which is also a Sturm–Liouville equation; in particular, $J_p(ay)$ satisfies this. Substituting in, we find $$ ( x J_p(\alpha x)')'J_p(\beta x) - (x J_p(\beta x)')' J_p(\alpha x) = \left( \frac{p^2}{x}-\alpha^2 x \right) J_p'(\alpha x))'J_p(\beta x) - \left( \frac{p^2}{x}-\beta^2 x \right) J_p(\alpha x) = (\beta^2-\alpha^2)x J_p(\alpha x)J_p(\beta x), $$ as required.
Exactly the same approach works for any eigenfunctions of a Sturm–Liouville equation: if $u_1$ and $u_2$ satisfy $-(pu_1')'+qu_1 = \lambda_1 u_1$ and $-(pu_2')'+qu_2 = \lambda_2 u_2$, then $$ \int wu_1u_2 = \frac{p (u_2'u_1-u_1'u_2)}{\lambda_1-\lambda_2} + C, \tag{$\dagger$} $$ since the derivative of the right-hand side's numerator is $$ (pu_2')'u_1+pu_2'u_1'-(pu_1')'u_2 - pu_1'u_2' = (q-\lambda_2 w)u_2u_1 - (q-\lambda_1 w)u_1u_2 = (\lambda_1-\lambda_2)wu_1u_2. $$
The slightly awkward thing here is that Bessel functions can be treated as solutions to the Sturm–Liouville equation (*) with two different eigenvalues: either $1/x$ or $x$ can be treated as the weight, with $p^2$ or $a^2$ being the eigenvalue. The latter is solved by ${\scr I}_p(ax)$ for fixed $p$ (where ${\scr I}$ denotes any solution to (*)), with $a$ taking values in the discrete set chosen so that the boundary conditions work), while the former is solved by varying $p$ and keeping $a$ fixed. The latter is also unusual in that $a$ is accounted for by simply scaling the argument, which is why the question's expression looks a bit different from the general version ($\dagger$) given above.