Motivation
I have the following non-homogeneous Bessel differential equation $$\frac{\mathrm{d}^2R}{\mathrm{d}r^2}+\frac1r\frac{\mathrm{d}R}{\mathrm{d}r}+\alpha^2R=J_0(\alpha r)$$
I want to find the general solution for this ODE. I know that the general solution can be written as
$$R=R_h+R_p$$
where $R_h$ is the basis for the homogeneous ODE and is known to be
$$R_h=C_1 \, J_0(\alpha r) + C_2 \, Y_0(\alpha r).$$
Next, for finding the particular solution $R_p$ one can use the method of variation of parameters to obtain a particular solution using the homogeneous ones. However, this will lead to some messy solution with hard integrals! As the MAPLE or WOLFRAM both use this technique their result is not valuable for me. I know that
$$R_p=\frac1{2 \alpha^2}\left[\alpha rJ_1(\alpha r)\right] \tag{*}$$
and it can be verified by putting it into the ODE. In fact, I saw this in some published paper.
Question
Is there an elegant method to obtain such a particular solution mentioned in $(*)$?