I am trying to learn a technique that an author uses for a boundary value problem and, so, I need to be able to use Mathematica to be able to solve the boundary conditions and then plot the final solution so that I can check alongside the journal article. The problem is that the equations involve multiple Bessel functions and the necessity to integrate them. For example, there are four unknown functions similar to the form:
$G(r) = \sum_{n=1}^{\infty} A_{i}h_{j}(\mu_{i}, r) + A_{-i}h_{j}(-\mu_{i}, r)$
where the $h_{j}$ will be linear combinations of the zero order Bessel function and the first order Bessel function and there are a total of 6 different $h_{j}$'s. Then, you have to take a vector of these functions and take the scalar product of it with a vector of the $h_{j}$'s and then finally integrate it over a radius to be able to find the unknown coefficients $A_{i}$ and $A_{-i}$. I attached a couple of screen shots of the relevant pages of one of the journal articles. However, my biggest problem in trying to duplicate the author's results are in using Mathematica to do the final integrals for the coefficients. If I let n=1, the functions will integrate nicely and I can then solve for the coefficients, but for n = 2 or greater, Mathematica stops integrating the functions nicely and just doesn't integrate at all. Then, if I try to NIntegrate, it gives me an error due to the coefficients being unknown. So, I was mostly looking to see if anyone had any advice in potentially making Mathematica able to integrate the equations fully, or, maybe somehow have NIntegrate factor out the coefficients from their grouped terms and then numerically integrate the result, so that I can still get a system of coefficient equations to then solve? I can't find anything related to the latter possibility and I figured that if Mathematica isn't solving the integrals analytically, then it may not be possible to find them analytically, forcing me to need to use NIntegrate. Any help would be appreciated.
g11[x_, r_] = -BesselJ[1, x * r];
g21[x_, r_] = BesselJ[0, x*r];
g31[x_, r_] = (c12 - c11) *x* r*BesselJ[1, x*r];
g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];
g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] -
x * r*BesselJ[0, x*r];
g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];
g32[x_, r_] = (c11 - c12) *x*
r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] -
x * r*BesselJ[0, x*r]);
g42[x_, r_] = (c11 - c12) *x*
r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] -
x* r*BesselJ[1, x*r]);
U1[x_] = -g31[x, 1]/g32[x, 1];
h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];
h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];
h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];
h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];
g11[x_, r_] = -BesselJ[1, x * r];
g21[x_, r_] = BesselJ[0, x*r];
g31[x_, r_] = (c12 - c11) *x* r*BesselJ[1, x*r];
g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];
g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] -
x * r*BesselJ[0, x*r];
g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];
g32[x_, r_] = (c11 - c12) *x*
r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] -
x * r*BesselJ[0, x*r]);
g42[x_, r_] = (c11 - c12) *x*
r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] -
x* r*BesselJ[1, x*r]);
U1[x_] = -g31[x, 1]/g32[x, 1];
h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];
h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];
h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];
h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];
Orz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m =
1\), \(1\)]\((Q[m]*h3[u[\([m]\)], r]\ + \
R[m]*h3[\(-u[\([m]\)]\), r])\)\)
Ozz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m =
1\), \(1\)]\((Q[m]*h4[u[\([m]\)], r]\ + \
R[m]*h4[\(-u[\([m]\)]\), r])\)\) + B0*\[Alpha]*r;
ur[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(1\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([m]\)]*l\)]*h1[u[\([m]\)], r]\ + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([m]\)]\)*l\)]*
h1[\(-u[\([m]\)]\), r])\)\) + B0*\[Beta]*r;
uz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(1\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([m]\)]*l\)]*h2[u[\([m]\)], r]\ + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([m]\)]\)*l\)]*
h2[\(-u[\([m]\)]\), r])\)\) + A0 + B0*l;
F0[r_] = {{0}, {0}, {Orz[r]}, {Ozz[r]}};
F1[r_] = {{ur[r]}, {uz[r]}, {0}, {P}};
Psi1T[m_, r_] = {{h1[x, r], h2[x, r], h3[x, r], h4[x, r]}};
Psi2T[m_, r_] = {{h1[-x, r], h2[-x, r], h3[-x, r], h4[-x, r]}};
Psi0T = {{0, 1, 0, 0}};
Psi01T[r_] = {{\[Beta]*r, 0, 0, \[Alpha]*r}};
(J = {{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0,
0}}) // MatrixForm
YY[x_, r_] = Psi2T[x, r].J.(F0[r] - E^(-x*l)*F1[r]);
XX[x_, r_] = Psi1T[x, r].J.(F0[r] - E^(x*l)*F1[r]);
DD[r_] = Psi0T.J.(F0[r] - F1[r]);
UUU[r_] = Psi01T[r].J.(F0[r] - F1[r]);
VVV[r_] = Psi0T.J.F0[r];
EE[r_] = UUU[r] - l*VVV[r];
P = 10;
l = 5;
Coeff1 = Integrate[YY[-u[[1]], r], {r, 0, 1}]
