Answer
(a) Using symmetry of the sphere, the potential depends only on the distance $r$ from $P$ to the center of the sphere. For convenience, we choose $P = \left( {0,0,r} \right)$ on the $z$-axis (with $r \ne R$) to compute $V\left( P \right)$.
(b) We show that
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi }}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{R^2} + {r^2} - 2Rr\cos \phi } }}$
(c) We show that
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left( {\left| {R + r} \right| - \left| {R - r} \right|} \right)$
(d) We verify Eq. (12) for $V$
Work Step by Step
(a) The gravitational potential due to $S$ at a point $P = \left( {a,b,c} \right)$ is equal to
$V\left( P \right) = - G\mathop \smallint \limits_{}^{} \mathop \smallint \limits_S^{} \dfrac{{\delta {\rm{d}}S}}{{\sqrt {{{\left( {x - a} \right)}^2} + {{\left( {y - b} \right)}^2} + {{\left( {z - c} \right)}^2}} }}$
Since the mass density is $\delta = \dfrac{m}{{4\pi {R^2}}}$, so
$V\left( P \right) = - \dfrac{{mG}}{{4\pi {R^2}}}\mathop \smallint \limits_{}^{} \mathop \smallint \limits_S^{} \dfrac{{{\rm{d}}S}}{{\sqrt {{{\left( {x - a} \right)}^2} + {{\left( {y - b} \right)}^2} + {{\left( {z - c} \right)}^2}} }}$
By the symmetrical structure of the sphere, the total contribution from each point $\left( {x,y,z} \right)$ on the surface of the sphere on $P$ can be considered as if the potential is due to the point in the center of the sphere on $P$. Therefore, the potential depends only on the distance $r$ from $P$ to the center of the sphere. For convenience, we can choose $P = \left( {0,0,r} \right)$ on the $z$-axis (with $r \ne R$) to compute $V\left( P \right)$.
(b) From part (a) we get
$V\left( P \right) = - \dfrac{{mG}}{{4\pi {R^2}}}\mathop \smallint \limits_{}^{} \mathop \smallint \limits_S^{} \dfrac{{{\rm{d}}S}}{{\sqrt {{{\left( {x - a} \right)}^2} + {{\left( {y - b} \right)}^2} + {{\left( {z - c} \right)}^2}} }}$
Since we choose $P = \left( {0,0,r} \right)$, that is $a=0$, $b=0$, $c=r$, the integral becomes
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi {R^2}}}\mathop \smallint \limits_{}^{} \mathop \smallint \limits_S^{} \dfrac{{{\rm{d}}S}}{{\sqrt {{x^2} + {y^2} + {{\left( {z - r} \right)}^2}} }}$
We express this integral using spherical coordinates, where
$x = R\cos \theta \sin \phi $, ${\ \ \ }$ $y = R\sin \theta \sin \phi $, ${\ \ \ }$ $z = R\cos \phi $
for $0 \le \theta \le 2\pi $, $0 \le \phi \le \pi $. In spherical coordinates, we have $dS = {R^2}\sin \phi d\theta d\phi $. So,
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi {R^2}}}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{{R^2}\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{{\left( {R\cos \theta \sin \phi } \right)}^2} + {{\left( {R\sin \theta \sin \phi } \right)}^2} + {{\left( {R\cos \phi - r} \right)}^2}} }}$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi }}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{R^2}{{\cos }^2}\theta {{\sin }^2}\phi + {R^2}{{\sin }^2}\theta {{\sin }^2}\phi + {R^2}{{\cos }^2}\phi - 2Rr\cos \phi + {r^2}} }}$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi }}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{R^2}{{\sin }^2}\phi + {R^2}{{\cos }^2}\phi - 2Rr\cos \phi + {r^2}} }}$
Hence,
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi }}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{R^2} + {r^2} - 2Rr\cos \phi } }}$
(c) From part (b) we get
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4\pi }}\mathop \smallint \limits_{\phi = 0}^\pi \mathop \smallint \limits_{\theta = 0}^{2\pi } \dfrac{{\sin \phi {\rm{d}}\theta {\rm{d}}\phi }}{{\sqrt {{R^2} + {r^2} - 2Rr\cos \phi } }}$
Write $u = {R^2} + {r^2} - 2Rr\cos \phi $. So, $du = 2Rr\sin \phi d\phi $. The integral becomes
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{8\pi Rr}}\mathop \smallint \limits_{\theta = 0}^{2\pi } \mathop \smallint \limits_{u = {R^2} + {r^2} - 2Rr}^{{R^2} + {r^2} + 2Rr} \dfrac{1}{{\sqrt u }}{\rm{d}}\theta {\rm{d}}u$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{8\pi Rr}}\mathop \smallint \limits_{\theta = 0}^{2\pi } {\rm{d}}\theta \mathop \smallint \limits_{u = {{\left( {R - r} \right)}^2}}^{{{\left( {R + r} \right)}^2}} \dfrac{1}{{\sqrt u }}{\rm{d}}u$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{4Rr}}\mathop \smallint \limits_{u = {{\left( {R - r} \right)}^2}}^{{{\left( {R + r} \right)}^2}} {u^{ - 1/2}}{\rm{d}}u$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left( {{u^{1/2}}|_{{{\left( {R - r} \right)}^2}}^{{{\left( {R + r} \right)}^2}}} \right)$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left[ {{{\left( {{{\left( {R + r} \right)}^2}} \right)}^{1/2}} - {{\left( {{{\left( {R - r} \right)}^2}} \right)}^{1/2}}} \right]$
Hence,
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left( {\left| {R + r} \right| - \left| {R - r} \right|} \right)$
(d) From part (c) we get
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left( {\left| {R + r} \right| - \left| {R - r} \right|} \right)$
1. Case $r>R$
$V\left( {0,0,r} \right) = - \dfrac{{mG}}{{2Rr}}\left( {R + r - \left( {r - R} \right)} \right) = - \dfrac{{mG}}{r}$
2. Case $r
