<s>[INST] Definite solution for the mean distance from an external point to the surface of a sphere. 
Sphere, radius $E$, is centred at point $O$ $[0,0,0]$.
External point $Q$ is at $[D,0,0]$. 
I can slice the sphere by making multiple planar cuts parallel to the $YZ$ plane to produce circular zones (quasi-discs) of equal infinitessimal width $dx$ and position $x_i$.  The curved surface of each zone has the same surface area.  The radius $R_i$ of any zone $i$ lies in the $YZ$ plane and has magnitude: $R_i = \sqrt{E^2 - x_i^2}  .$
The distance $L_i$ from point $Q$ to any point $P_i$ in the zone $i$ is given by:-
$$L_i = \sqrt{ (D-x_i)^2 + R_i^2} 
 = \sqrt{ D^2 -2Dx_i + x_i^2 + E^2 - x_i^2} 
= \sqrt{ D^2 -2Dx_i + E^2 } $$
$$L_i  = D \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} }. $$
Now I wish to integrate this expression over the range $-E\cdots+E$ and then divide by $2E$ to obtain the average value $\bar{L}$ thus
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx.$$
I can make a Taylor Series approximation of the square root term using the standard formula $(1+x)^{0.5} = 1 + x/2 - x^2/8 + x^3/16 -\cdots$  
Applying this to the square root term I obtain
$$
\sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } = 1-x_i/D + (1/2)E^2/D^2 -(1/2)x_i^2/D^2 + (1/2)x_iE^2/D^3  + \cdots$$ where the subsequent terms on the RHS diminish in magnitude.  By dropping terms with odd powers of $x_i$ (because
they will go to zero when integrating over the range $-E \le x_i \le +E$ ) and consolidating and integrating I come up with the following approximate result:-
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx
\approx
D\left(1 + \frac{E^2}{3D^2} + A\right) $$
where A is a small term whose expression depends on the number of terms evaluated in the Taylor series approximation.  From numerical modelling it appears plausible that the term $A$ vanishes if the Taylor Series is extended to infinite terms and thus the definite solution would be given by:-
$$\bar{L} = 
D\left(1 + \frac{E^2}{3D^2} \right)
= D + \frac{E^2}{3D}
$$
MY QUESTION
Is there a way to derive a definite solution to this problem?
 [/INST] Following suggestions from coffeemath & user2566092 I found the following standard rule:-
$$
\int (ax + b) ^{p/n} = \frac{n}{(n+p)a} (ax + b) ^{1+p/n} +C
 $$
for $ p = \pm1, \pm2, \cdots p \ne -n .$
(Source: The Universal Encyclopedia of Mathematics, page 590; a similar equation is here).
Applying that to the expression for $\bar {L}$
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx
=
 \frac{D}{2E} \int_{-E}^{+E} \left[ \left(\frac{-2}{D}\right)x_i  +    \left(1 + \frac{E^2}{D^2}\right)  \right]^{1/2} dx$$
whence $a=-2/D,b = 1 +(E^2/D^2), p = 1, n = 2$ and so
$$\bar{L} = 
 \frac{D}{2E} 
\frac{2}{3(-2/D)}
\left[ 
\left(
\left(\frac{-2}{D}\right)x_i  +    \left(1 + \frac{E^2}{D^2}\right)  \right)^{3/2}+constant\right]_{-E}^{+E} $$
then cancelling and dropping the constant leads to
$$\bar{L} = 
\frac{-D^2}{6E}
\left[ \left(
\frac{-2Dx_i+D^2+E^2}{D^2}   \right)^{3/2}\right]_{-E}^{+E} 
 = 
\frac{-D^2}{6E} \left( \frac{1}{D^2}\right)
^{3/2}
\left[ \left(
-2Dx_i+D^2+E^2  \right)^{3/2}\right]_{-E}^{+E} $$
expanding and shuffling minus signs
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
\left(2DE+D^2+E^2  \right)^{3/2}
-
\left(-2DE+D^2+E^2  \right)^{3/2}
\right]
$$
this simplifies to
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
\left( (D+E)^2  \right)^{3/2}
-
\left((D-E)^2 \right)^{3/2}
\right]
 = 
\frac{1}{6DE}
\left[ 
 (D+E)^3  -(D-E)^3 
\right]
$$
and then
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
 6D^2E +2E^3 
\right]
=
D + \frac{E^2}{3D}
$$
and finally
$$\bar{L} = 
D \left( 1+ \frac{E^2}{3D^2} \right).
$$
</s>