Calculate the field of a sphere with a homogenous distributed charge $\rho_0$ using Gaussian theorem and the 1st Maxwell Equation:
$$ \oint \limits_{S} \vec{E} \cdot d \vec{S} \underset{\substack{\uparrow \\ \text { Gauss }}}{=} \int \limits_{V} (\vec{\nabla} \cdot \vec{E}) dV \underset{\substack{\uparrow \\ \text { 1. Maxwell }}}{=} \frac{1}{\epsilon_{0}} \int \limits_{V} \rho dV=\frac{Q(V)}{\varepsilon_{0}} $$
I don't have problems with the volume integral:
Let $R$ be the radius of the sphere:
$$Q(V)=\int \limits_{V} \rho dV=\rho_o \int \limits_{0}^{R} r^2 d r\int \limits_{0}^{2 \pi} d \phi \int \limits_{0}^{\pi} d \theta \sin \theta= \rho_0 \frac{4\pi }{3}R^3 $$
In spherical coordinates the field is given by $\vec{E}(\vec{r}) =E_{r} \vec{e}_{r}+E_{\theta} \vec{e}_{\theta}+E_{\phi} \vec{e}_{\phi}$
I thing because of the radial symmetry, I can assume: $$\vec{E}(\vec{r}) =E_{r} \vec{e}_{r}$$
with $\vec{e}_{r}=[(\sin \theta) (\cos \phi),(\sin \theta) (\sin \phi),\cos \theta]^t $.
Our convention is that the normal surface vector points outward i.e. $d \vec{S}\| \vec{e}_{r}$.
How do I calculate: $\oint \limits_{S} \vec{E} \cdot d \vec{S}$?
Let's assume the sphere is centered at $0$ and has radius $R$. Then the outward pointing normal vector at $\vec{x}\in S$ is $\vec{x}/R$ (since $\vec{x}$ has lenght $R$). You want to find $\vec{E}$, which, as you argued by spherical symmetry, has constant magnitude on the entire sphere. Thus $$\int_S \vec{E}\cdot d\vec{S} = \int_S \vec{E}\cdot\frac{\vec{x}}{R} dS = \int_S E dS =E\int_SdS=E\cdot 4πR^2,$$ where, in your notation, $E = E_r$.