While solving a static problem regarding an isotropic cylinder under its own weight, I ran into the following integral:
$$ I_0 = \int_0^a r^2 J_0(\mu_ir) J_1(\mu_jr)dr \tag 1$$
where $r$ represents the radial coordinate, $J_0$ and $J_1$ are the Bessel functions of the first kind with order zero and one, $a$ is the cylinder's radius ($a = r_{max}$) and $\mu_n$ is the $n$-th root of the following equation:
$$\Bigg[ (\mu a)^2 - { 2 c_{11} \over c_{11} + c_{12} } \Bigg]J_1^2(\mu a) + (\mu a)^2 J_0^2(\mu a) = 0 \tag 2$$
In the above expression, $c_{11}$ is the longitudinal modulus and $c_{12}$ is the Lamé's constant. Some specific values of $\mu a$ that satisfy equation $(2)$ are given in the picture below.
I tried to solve the integral $(1)$, but it appears I went wrong somewhere. I am hoping somebody here can point me in the right direction so I can find my mistake. Here is what I did.
Step $\textbf 1$ - Solve the special case where $\mu_i = \mu_j = \mu$
The solution in this case is:
$$ \int_0^a r^2 J_0(\mu r) J_1(\mu r)dr = {a^2 \over 2 \mu} J_1^2(\mu a) \tag 3$$
which I verified using Wolfram Alpha.
Step $\textbf 2$ - Solve the case where $\mu_i \neq \mu_j$
Substep $1$ - Use the substitution method to reformulate equation $(1)$
Let $\mu_j r = x$, $\mu_i / \mu_j = k$ and $\mu_j a = b$. Now equation $(1)$ can be rewritten in the form:
$$ I_0 = {1 \over \mu_j ^3}\int_0^b x^2 J_0(kx) J_1(x)dx \tag 4$$
Now I need to solve the integral:
$$ I = \int_0^b x^2 J_0(kx) J_1(x)dx \tag 5$$
Substep $2$ - Use the Feynman techinique to solve equation $(5)$
Assuming that $x^2 J_0(kx) J_1(x)$ is continuous for all values of $k$ in the bounds of $[0, b]$, I used the Leibniz rule to swap the order of differentiation and integration:
$${dI(k) \over dk} = {d \over dk}\int_0^b x^2 J_0(kx) J_1(x)dx \tag 6$$ $${dI(k) \over dk} = \int_0^b x^2 {\partial J_0(kx) \over \partial k} J_1(x)dx \tag 7$$
$${dI(k) \over dk} = \int_0^b -x^3 J_1(kx) J_1(x)dx \tag 8$$
Note that I used the Bessel function properties to derive that $\partial J_0(kx) / \partial k$ is equal to $-x J_1(kx)$.
Substep $3$ - Use integration by parts to solve equation $(8)$
Let $u = -x^2$ and $dv = xJ_1(kx)J_1(x)dx$. It is easy to derive $du = -2x dx$. However, to determine $v$, I used Wolfram Alpha:
$$v = {x J_0(x)J_1(kx) - kx J_1(x) J_0(kx) \over k^2 - 1} \tag 9$$
Now, equation $(8)$ can be rewritten as:
$${dI(k) \over dk} = -x^2{x J_0(x)J_1(kx) - kx J_1(x) J_0(kx) \over k^2 - 1} \Bigg |_0^b - {\int_0^b} {x J_0(x)J_1(kx) - kx J_1(x) J_0(kx) \over k^2 - 1}(-2)x dx \tag {10} $$
Rearranging the terms in the above equation, and using the expression $(5)$ yields:
$${dI(k) \over dk} = {b^3J_0(b)J_1(kb) - b^3kJ_0(kb)J_1(b)\over 1-k^2} + {2k \over 1-k^2}I - {2 \over 1-k^2} {\int_0^b} x^2 J_0(x) J_1(kx) dx \tag {11}$$
Substep $4$ - Use integration by parts to reformulate $(11)$ in terms of $I$
Let $u = xJ_1(kx)$ and $dv = xJ_0(x)dx$. Using Bessel function properties yields $du = kxJ_0(kx)dx$ and $v = xJ_1(x)$. Now we can write:
$${\int_0^b} x^2 J_0(x) J_1(kx) dx = x^2 J_1(x) J_1(kx) \Bigg |_0^b - k \int_0^b x^2 J_0(kx) J_1(x) dx \tag {12}$$
$$ {\int_0^b} x^2 J_0(x) J_1(kx) dx = b^2 J_1(b) J_1(kb) - kI \tag {13}$$
Inserting the above expression into equation $(11)$ gives us the following differential equation:
$${dI(k) \over dk} = {b^3J_0(b)J_1(kb) - b^3kJ_0(kb)J_1(b)\over 1-k^2} + {2k \over 1-k^2}I - {2 \over 1-k^2} \Big( b^2 J_1(b) J_1(kb) - kI \Big) \tag {14}$$
which can be rewritten as:
$$ {dI(k) \over dk} - {4k \over 1 - k^2}I = -b^2 \Big( 2J_1(b) - bJ_0(b) \Big) {J_1(kb) \over 1 - k^2} - b^3 J_1(b) {k J_0(kb) \over 1 - k^2}\tag {15}$$
The above expresion can be simplified using the property $bJ_2(b) = 2J_1(b) - bJ_0(b)$:
$$ {dI(k) \over dk} - {4k \over 1 - k^2}I = -b^3 J_2(b) {J_1(bk) \over 1 - k^2} - b^3 J_1(b) {k J_0(bk) \over 1 - k^2} \tag {16}$$
Substep $5$ - Use the integration factor to solve the differential equation $(16)$
Substituting $I(k) = (1-k^2)^{-2}f(k)$ into the expression $(16)$ yields the following separable differential equation:
$${df(k) \over dk} = -b^3 J_2(b)(1 - k^2)J_1(bk) - b^3 J_1(b)(1 - k^2)kJ_0(bk) \tag {17}$$
$$\int df= \int \Bigg( -b^3 J_2(b)(1 - k^2)J_1(bk) - b^3 J_1(b)(1 - k^2)kJ_0(bk) \Bigg) dk \tag {18}$$
$$ f = -b^3 J_2(b) \int J_1(bk) dk +b^3 J_2(b) \int k^2 J_1(bk) dk - b^3 J_1(b) \int k J_0(bk) dk + b^3 J_1(b) \int k^3 J_0(bk) dk \tag {19}$$
Neglecting the integration constants, the solutions to the first, second, third and fourth integral are:
$$\int J_1(bk) dk = {1 - J_0(bk) \over b} \tag {20}$$ $$\int k^2 J_1(bk) dk = {k^2 J_2(bk) \over b} \tag {21}$$ $$\int k J_0(bk) dk = {kJ_1(bk) \over b} \tag {22}$$ $$\int k^3 J_0(bk) dk = {2k^2 J_2(bk) - b k^3 J_3(bk)\over b^2}\tag {23}$$
which means:
$$ f = -b^3 J_2(b) {1 - J_0(bk) \over b} + b^3 J_2(b) {k^2 J_2(bk) \over b} - b^3 J_1(b) {kJ_1(bk) \over b} + b^3 J_1(b) {2k^2 J_2(bk) - b k^3 J_3(bk)\over b^2} + C \tag {24}$$
where $C$ is the integration constant. Inserting the above expression into $I(k) = (1-k^2)^{-2}f(k)$, and rearranging the terms gives us the following function:
$$I = {1 \over (1-k^2)^2} \Bigg( C - b^3 J_2(b){1 - J_0(bk) - k^2 J_2(bk) \over b} - b^3 J_1(b) { bkJ_1(bk)-2k^2J_2(bk)+bk^3J_3(bk) \over b^2} \Bigg) \tag {25}$$
Substep $6$ - Determine the value of the constant $C$ by solving the equation $(5)$ for the value of $k = 0$
The value of $I(0)$ is:
$$ I(0) = \int_0^b x^2 J_1(x)dx = b^2 J_2(b) \tag {26}$$
Inserting $k=0$ and $I(0) = b^2 J_2(b)$ into equation $(25)$ gives us the value of the integration constant as:
$$ C = b^2 J_2(b) \tag {27}$$
which gives us the function:
$$I = {1 \over (1-k^2)^2} \Bigg( b^2 J_2(b) - b^3 J_2(b){1 - J_0(bk) - k^2 J_2(bk) \over b} - b^3 J_1(b) { bkJ_1(bk)-2k^2J_2(bk)-bk^3J_3(bk) \over b^2} \Bigg) \tag {28}$$
Rearraning the terms in the above equation gives us:
$$ I = {b \over (1 - k^2)^2} \Bigg[ bJ_2(b) \Big( J_0(bk) + k^2 J_2(bk) \Big) + J_1(b) \Big( 2k^2J_2(bk)-bkJ_1(bk)-bk^3J_3(bk) \Big) \Bigg] \tag{29} $$
Substep $7$ - Do the resubstitution and determine the final expression of $I_0$
Since $I_0 = {1 \over \mu_j ^3}I(\mu_i/\mu_j)$, we can write it down as:
$$ I_0 = {a \over \mu_j^2\Big(1 - {\mu_i^2 \over \mu_j^2}\Big)^2} \Bigg[ \mu_jaJ_2(\mu_ja) \Big( J_0(\mu_ia) + {\mu_i^2 \over \mu_j^2} J_2(\mu_ia) \Big) + J_1(\mu_ja) \Big( 2{\mu_i^2 \over \mu_j^2}J_2(\mu_ia)-\mu_iaJ_1(\mu_ia)-\mu_ja {\mu_i^3 \over \mu_j^3}J_3(\mu_ia) \Big) \Bigg] \tag{30} $$
Now, we get the final expression as:
$$ I_0 = { a \over \Big(\mu_j^2 - {\mu_i^2}\Big)^2} \Bigg[ \mu_jaJ_2(\mu_ja) \Big( \mu_j^2J_0(\mu_ia) + {\mu_i^2} J_2(\mu_ia) \Big) + \mu_iJ_1(\mu_ja) \Big( 2{\mu_i}J_2(\mu_ia)-\mu_j^2aJ_1(\mu_ia)-{\mu_i^2}aJ_3(\mu_ia) \Big) \Bigg] \tag{31} $$
I tried to verify my results by evaluating numerically the integral $I_0$ using the midpoint rule in Matlab. For the case of $\mu_i = \mu_j=\mu = 4$ and $a=2$, the analitical solution $(3)$ had the same value as my numerical evaluation (which is $I_0 = 0.02753$). But for the case where $\mu_i = 4$, $\mu_j = 5$ and $a=2$, the results did not match. The numerical value I evaluated was $I_0 = -0.10535$, while the analytical was $I_0 = 0.14188$. That means I made a mistake somewhere, but I don't know where. Note that I am aware that I used real values instead of complex values for $\mu_i$ and $\mu_j$. That is just for numerical verification.
It would be of help to me if someone pointed me in the right direction to solve this problem.
EDIT: I apolagise, the error was in my Matlab code. The solution I derived seems correct. I will close this question.