I'd like to start with the original problem statement and work from there, since there are a few small mistakes in that integral you've got up there. This is Problem 2.7 in Griffiths' Introduction to Electrodynamics (Fourth Edition).
Let $\mathbf{r}$ be the vector from any point on the sphere's surface to $P$ (because that script $\mathbf{r}$ that Griffiths uses is annoyingly hard to type in MathJax), and let $r$ be its length. Also, let $\alpha$ be the angle between $\mathbf{r}$ and the $z$-axis, and let $\phi$ be the azimuthal angle. Then, as in the figure, $\theta$ is the polar angle, $R$ is the radius of the sphere, and $z$ is the distance from the center of the sphere to $P$.
Since the surface of the sphere has a total charge $q$ distributed with a uniform density $\sigma$, each infinitesimal part $dA$ of the sphere's surface carries an infinitesimal piece of charge $dq=\sigma dA$. Expressing that in spherical coordinates, we have $dq=\sigma R^2 \sin\theta\ d\theta\ d\phi$.
We can see from the symmetry of the problem setup that any horizontal components of the electric field at $P$ will cancel out, so $\mathbf{E}=E_z\hat{\mathbf{z}}$, which simplifies things a bit -- we only need to find $E_z$.
Now, Coulomb's law tells us that the strength of the electric field generated by each infinitesimal piece of charge $dq$ is
$$dE=\frac1{4\pi\epsilon_0}\frac{dq}{r^2}.$$
More specifically, in this case (at $P$),
$$dE=dE_z=\frac1{4\pi\epsilon_0}\frac{dq}{r^2}\cos\alpha,$$
where we need the $\cos\alpha$ term to represent the contribution to the electric field at $P$ from all points other than the one directly between $P$ and the center of the sphere.
Since we have $dE_z$ at the moment, the logical thing to do is to integrate to find $E_z$, but we're going to need to re-express some terms first. $dq$ (as we've already seen), $r^2$, and $\cos\alpha$ can all be re-expressed in terms of combinations of variables over which we're going to integrate and constants.
By the law of cosines,
$$r^2=R^2+z^2-2Rz\cos\theta.$$
We can also re-express $\cos\alpha$ in a similar way:
$$\begin{align}
R^2&=r^2+z^2-2rz\cos\alpha \\
\cos\alpha&=\frac{R^2-r^2-z^2}{-2rz} \\
&=\frac{r^2+z^2-R^2}{2rz} \\
&=\frac{R^2+z^2-2Rz\cos\theta+z^2-R^2}{2z\sqrt{R^2+z^2-2Rz\cos\theta}} \\
&=\frac{2z(z-R\cos\theta)}{2z\sqrt{R^2+z^2-2Rz\cos\theta}} \\
&=\frac{z-R\cos\theta}{\sqrt{R^2+z^2-2Rz\cos\theta}}
\end{align}$$
Pulling everything together, we have:
$$\begin{align}
dE_z&=\frac1{4\pi\epsilon_0}\frac{\sigma R^2 \sin\theta\ d\theta\ d\phi}{R^2+z^2-2Rz\cos\theta}\frac{z-R\cos\theta}{\sqrt{R^2+z^2-2Rz\cos\theta}} \\
&=\frac1{4\pi\epsilon_0}\frac{\sigma R^2 \sin\theta(z-R\cos\theta)}{(R^2+z^2-2Rz\cos\theta)^{3/2}}d\theta\ d\phi
\end{align}$$
Now we can integrate. Since $\theta$ is the polar angle, we have $0\le\theta\le\pi$, and since $\phi$ is the azimuthal angle, we have $0\le\phi\le2\pi$; note that the definitions of $\theta$ and $\phi$ I use are the opposite of those you used in your integral. Note also that since $dr$ never appeared in our earlier calculations, there's no need to integrate over $r$.
Rearranging a little bit, our integral emerges at last:
$$E_z=\frac1{4\pi\epsilon_0}\sigma R^2\int_0^{2\pi}\int_0^{\pi}\frac{\sin\theta(z-R\cos\theta)}{(R^2+z^2-2Rz\cos\theta)^{3/2}}d\theta\ d\phi$$
That's a mess. Let's see what we can do about that.
It's easy to integrate over $\phi$ first, since $\phi$ doesn't appear anywhere in the integrand; that simplifies the integral a little bit, to this:
$$E_z=\frac1{2\epsilon_0}\sigma R^2\int_0^{\pi}\frac{\sin\theta(z-R\cos\theta)}{(R^2+z^2-2Rz\cos\theta)^{3/2}}d\theta$$
We can now make the substitution $u=z-R\cos\theta$, and once we do that, a few interesting things happen:
$$\begin{align}
du&=R\sin\theta\ d\theta \implies \frac1Rdu=\sin\theta\ d\theta \\
\theta&=0 \implies u=z-R\cos(0)=z-R \\
\theta&=\pi \implies u=z-R\cos(\pi)=z+R \\
u&=z-R\cos\theta \implies R\cos\theta=z-u
\end{align}$$
Plugging all of that back into the integral, we have:
$$\begin{align}
E_z&=\frac1{2\epsilon_0}\sigma R^2\int_{z-R}^{z+R}\frac{\frac1Ru}{(R^2+z^2-2z(z-u))^{3/2}}du \\
&=\frac1{2\epsilon_0}\sigma R\int_{z-R}^{z+R}\frac{u}{(R^2+z^2-2z^2+2zu)^{3/2}}du \\
&=\frac1{2\epsilon_0}\sigma R\int_{z-R}^{z+R}\frac{u}{(R^2-z^2+2zu)^{3/2}}du
\end{align}$$
Let's make another substitution now: $v=R^2-z^2+2zu$. This has similar effects to the last substitution:
$$\begin{align}
u&=\frac1{2z}(v-R^2+z^2) \\
dv&=2z\ du \implies \frac1{2z}dv=du
\end{align}$$
$$\begin{align}
u=z-R \implies v&=R^2-z^2+2z(z-R) \\
&=R^2-z^2+2z^2-2zR \\
&=R^2+z^2-2zR \\
&=(R-z)^2 \\
u=z+R \implies v&=R^2-z^2+2z(z+R) \\
&=R^2-z^2+2z^2+2zR \\
&=R^2+z^2+2zR \\
&=(R+z)^2
\end{align}$$
Our integral now looks like this:
$$\begin{align}
E_z&=\frac1{2\epsilon_0}\sigma R\int_{(R-z)^2}^{(R+z)^2}\frac{\frac1{2z}(v-R^2+z^2)}{v^{3/2}}\frac1{2z}dv \\
&=\frac1{8z^2\epsilon_0}\sigma R\int_{(R-z)^2}^{(R+z)^2}\frac{v-R^2+z^2}{v^{3/2}}dv \\
&=\frac1{8z^2\epsilon_0}\sigma R\int_{(R-z)^2}^{(R+z)^2}\left[v^{-1/2}-(R^2-z^2)v^{-3/2}\right]dv
\end{align}$$
Much more manageable, right? At long last, we can evaluate it!
$$\begin{align}
E_z&=\frac1{8z^2\epsilon_0}\sigma R\left[2v^{1/2}+2(R^2-z^2)v^{-1/2}\right]\Bigg|_{(R-z)^2}^{(R+z)^2} \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[v^{1/2}+(R^2-z^2)v^{-1/2}\right]\Bigg|_{(R-z)^2}^{(R+z)^2} \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[\sqrt{(R+z)^2}+\frac{R^2-z^2}{\sqrt{(R+z)^2}}-\sqrt{(R-z)^2}-\frac{R^2-z^2}{\sqrt{(R-z)^2}}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[R+z+\frac{R^2-z^2}{R+z}-\sqrt{(R-z)^2}-\frac{R^2-z^2}{\sqrt{(R-z)^2}}\right]
\end{align}$$
I agree that it looks a little weird to leave the $(R-z)^2$ terms under the square roots, but remember the second part of the hint from the problem statement? That's going to come into play now. We have two distinct cases, where $z\lt R$ (in other words, where $P$ is inside the sphere), and where $z\gt R$ (where $P$ is outside the sphere). The value of $\sqrt{(R-z)^2}$ will be different in each case, since we always want to avoid taking the square root of a negative number.
If $P$ is inside the sphere $(z\lt R)$:
$$\begin{align}
E_z&=\frac1{4z^2\epsilon_0}\sigma R\left[R+z+\frac{R^2-z^2}{R+z}-\sqrt{(R-z)^2}-\frac{R^2-z^2}{\sqrt{(R-z)^2}}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[R+z+\frac{R^2-z^2}{R+z}-R+z-\frac{R^2-z^2}{R-z}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2z+\frac{(R^2-z^2)(R-z)-(R^2-z^2)(R+z)}{R^2-z^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2z+\frac{(R^2-z^2)(R-z-R-z)}{R^2-z^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2z-\frac{2z(R^2-z^2)}{R^2-z^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left(2z-2z\right) \\
&=0
\end{align}$$
If $P$ is outside the sphere $(z\gt R)$:
$$\begin{align}
E_z&=\frac1{4z^2\epsilon_0}\sigma R\left[R+z+\frac{R^2-z^2}{R+z}-\sqrt{(R-z)^2}-\frac{R^2-z^2}{\sqrt{(R-z)^2}}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[R+z+\frac{R^2-z^2}{R+z}-z+R-\frac{R^2-z^2}{z-R}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2R+\frac{(R^2-z^2)(z-R)-(R^2-z^2)(R+z)}{z^2-R^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2R+\frac{(R^2-z^2)(z-R-R-z)}{z^2-R^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2R-\frac{2R(R^2-z^2)}{z^2-R^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left[2R+\frac{2R(z^2-R^2)}{z^2-R^2}\right] \\
&=\frac1{4z^2\epsilon_0}\sigma R\left(2R+2R\right) \\
&=\frac1{4z^2\epsilon_0}\sigma R\left(4R\right) \\
&=\frac1{z^2\epsilon_0}\sigma R^2 \\
&=\frac1{4\pi\epsilon_0}\frac{q}{z^2} \qquad\qquad (\text{since}\ q=4\pi R^2\sigma)
\end{align}$$
And there you have it!
$$\mathbf{E}=
\begin{cases}
\mathbf{0}, & z\lt R \\[1ex]
\displaystyle{\frac1{4\pi\epsilon_0}\frac{q}{z^2}}\hat{\mathbf{z}}, & z\gt R
\end{cases}$$
