AssumpCS={CS->L*(L+hd)}
AssumpCS2={CS2->CS-(hd+1)}
AssumpCS3={CS3->CS2-(hd+1)}
AssumpMs={Ms[r]->CS CS2 (hd^2)-2 hd (1+hd) CS2 (r^2) (ω^2)+(((1+hd))^2) (r^4) (ω^4)}
α=L+hd/2
Print["Bessel functions to next-to-next-to-leading-order"]
besselj[x_]=Simplify[1-((x^2)/(2*(2α+2)))+((x^4)/(4*2*(2α+2)(2α+4)))]
AssumpBesselEq={j''[ω r]->FullSimplify[-(((2 α+1)/(ω r))j'[ω r]+j[ω r])]}
AssumpBesselSub={j[ω r]->besselj[ω r],j'[ω r]->besselj'[ω r]}

Trr=(1/(Nld*Omd))∫ (Subscript[N, abL])[Om](Subscript[T, ab])[r]\[DifferentialD]Om
Print["pref:"]
pref=((Nld*Omd*(hd+1)*(r^(L+hd+3)))/(2 CS))
AssumpXYZ={X[r]->(((ω^2)(r^(hd+5)))/(Ms[r]))Trr,Y[r]->(((ω^2)(r^3)(CS CS2 (hd-1)(hd^2)- 2 hd ((( hd+1))^2)CS2 (ω^2)(r^2)+(hd+3)((( hd+1))^2)(ω^4)(r^4)))/(((Ms[r]))^2)) Trr, Z[r]->(1/(((Ms[r]))^2))((CS^(2 ))(CS2^2)(hd^3)
-CS CS2 (hd^2)(hd+1)(2CS+(hd+1)((hd^2)+hd-4))(ω^2)(r^2)+CS2 hd ((( hd+1))^2)(CS+2hd(hd+1)(hd+3))(ω^4)(r^4)-(((hd+1))^4)(hd+2)(hd+3)(ω^6)(r^6))Trr}
Print["J:"]
J=Z[r]-(((hd+1))^2)(∂_r)Y[r]-(((hd+1))^2)(∂_r)((1/(r^(hd+1)))(∂_r)X[r])
Print["𝒯:"]
𝒯=pref*((hd/r)(∂_r)J+(ω^2)J-((hd CS3)/((hd+1)(r^2)))J);
Print["original integrand 𝒯 j[ω r]:"]
integrand=FullSimplify[𝒯 j[ω r]]
fullderivative1=(∂_r)(pref * (hd/r) * j[ω r]*J);
fullderivative2=(∂_r)(pref* (r^(-(hd+4))) ((hd (CS3+(1+hd) (2+hd+L))-(1+hd) (r^2) (ω^2)) j[r ω]+hd (1+hd) r ω j'[r ω]) (1+hd) (-(1+hd) X[r]+r ((r^(1+hd)) Y[r]+X'[r])));
Print["After integrations by parts:"]
integrand2=FullSimplify[integrand-fullderivative1-fullderivative2]
Print["After Bessel Equation:"]
integrand3=FullSimplify[integrand2/.AssumpBesselEq]

fullderivative3=(∂_r)(pref*(r^(-(hd+5)))(1+hd)(-r (hd L (CS3+(1+hd) (2+hd+L))-(1+hd) (2+hd+L) (r^2) (ω^2)) j[r ω]+(r^2) ω (-hd (2+CS3+2 hd)+(1+hd) (r^2) (ω^2)) j'[r ω])X[r]);
Print["After 2nd integrations by parts:"]
integrand22=FullSimplify[integrand3-fullderivative3]
Print["After 2nd Bessel Equation:"]
integrand32=FullSimplify[integrand22/.AssumpBesselEq]

Print["After plugging X,Y,Z, without common prefactor:"]
commonfactor=(1/(2 CS )) (r^(1+hd+L))(∫ (Subscript[N, abL])[Om](Subscript[T, ab])[r]\[DifferentialD]Om);
Print["work on Trr[r] j'[r ω]:"]
Factorjderv=FullSimplify[((FullSimplify[FullSimplify[integrand32/.{j[ω r]->0}]/.AssumpXYZ])/commonfactor)/.AssumpCS3]
Factorjderv2=FullSimplify[FullSimplify[Factorjderv/.AssumpMs]/.AssumpCS2]
Print["Final Trr[r] j'[r ω]:"]
Factorjderv3=FullSimplify[Factorjderv2/.AssumpCS]
Print["work on Trr[r] j[r ω]:"]
Factorj=FullSimplify[((FullSimplify[FullSimplify[integrand32/.{j'[r ω]->0}]/.AssumpXYZ])/commonfactor)/.AssumpCS3]
Factorj2=FullSimplify[FullSimplify[Factorj/.AssumpMs]/.AssumpCS2]
Factorj3=FullSimplify[Factorj2/.AssumpCS]
Print["Final Trr[r] j[r ω]:"]
Factorj4=FullSimplify[Factorj3*((CS+(hd+1)(L+hd-(ω^2)(r^2)))/(hd+(hd^2)+L+2 hd L+(L^2)-(1+hd) (r^2) (ω^2)))]
Print["test:"]
FullSimplify[Factorj4-Factorj3/.AssumpCS]

Print["Final integrand:"]
integrand5=FullSimplify[commonfactor*(Factorjderv3+Factorj4)]
Print["test:"]
FullSimplify[FullSimplify[FullSimplify[FullSimplify[integrand32-integrand5/.AssumpXYZ]/.AssumpMs]/.AssumpCS2]/.AssumpCS]
integrand6=FullSimplify[integrand5/.{(Subscript[N, abL])[Om]->(((Subscript[x, abL])[Om,r])/(r^(L+2)))}]
Q=FullSimplify[∫(integrand6) \[DifferentialD]r]