R = 3.9
x = Cos[[Theta]]
b = (r^2 + R^2 - 2rRx)^0.5
NIntegrate[(Sqrt[15]/4)
Integrate[(R - rx)(3 x^2 - 1)E^(-1.5566b), {x, -1, 1}]r^4
E^(-3.313066*r), {r, 0, R}]
Your “inner” integrand:
(R - r*x)*(3 x^2 - 1)*E^(-1.5566*b)
evaluates to an expression that depends on theta and r. You probably want to integrate over theta. Something like (you will want to check this assumption so that it corresponds to your problem:
-Integrate[(R - r*x)*(3 x^2 - 1)*E^(-1.5566*b) D[x, \[Theta]], {\[Theta], 0, Pi}, Assumptions -> 0 < r < R]
This may not produce a closed form. So, let’s define a function:
integrand = (R - r*x)*(3 x^2 - 1)*E^(-1.5566*b) D[x, \[Theta]]
innerInt[rvar_?NumericQ] :=
With[{tmp = integrand /. r -> rvar}, -NIntegrate[ tmp, {\[Theta], 0, Pi}] ]
test it:
innerInt[.4]
it works.
Use it:
NIntegrate[(Sqrt[15]/4)*innerInt[r]*r^4*E^(-3.313066*r), {r, 0, R}]
This works too.
Again, you will want to check what I did against what you intended. But the steps should work.
If d Cos[u] is not accepted as integration-variable you should change it to ( - Sin[u] du ). You have to change the boundaries as well
Then: your inner integrand is really complicated containing transcedent functions.
So I tried it numerically in total:
f1[R_, r_?NumericQ] :=
NIntegrate[(R - r*Cos[u])*(3 Cos[u]^2 - 1)*
E^(-1.5566*(r^2 + R^2 - 2*r*R*Cos[u])^0.5) (-Sin[u]), {u, -Pi, 0}]
and
NIntegrate[(Sqrt[15]/4)*f1[3.9, r]*r^4*E^(-3.313066*r), {r, 0, 3.9}]
