How is Gauss' Law (integral form) arrived at from Coulomb's Law, and how is the differential form arrived at from that?

Let us for simplicity consider $n$ point charges $q_1$, $\ldots$, $q_n$, at positions $\vec{r}_1$, $\ldots$, $\vec{r}_n$, in the electrostatic limit, with vacuum permittivity $\epsilon_0$.

Now let us sketch one possible strategy to prove Gauss' law from Coulomb's law:

1. Deduce from Coulomb's law that the electric field at position $\vec{r}$ is $$\tag{1} \vec{E}(\vec{r})~=~ \sum_{i=1}^n\frac{q_i }{4\pi\epsilon_0}\frac{\vec{r}-\vec{r}_i}{|\vec{r}-\vec{r}_i|^3} . $$

2. Deduce the charge density $$\tag{2} \rho(\vec{r})~=~\sum_{i=1}^n q_i\delta^3(\vec{r}-\vec{r}_i). $$

3. Recall the following mathematical identity $$\tag{3}\vec{\nabla}\cdot \frac{\vec{r}}{|\vec{r}|^3}~=~4\pi\delta^3(\vec{r}) .$$ (This Phys.SE answer may be useful in proving eq.(3), which may also be written as $\nabla^2\frac{1}{|\vec{r}|}=-4\pi\delta^3(\vec{r})$).

4. Use eqs. (1)-(3) to prove Gauss' law in differential form $$\tag{4} \vec{\nabla}\cdot \vec{E}~=~\frac{\rho}{\epsilon_0} .$$

5. Deduce Gauss' law in integral form via the divergence theorem.

@Qmechanic's already provided a nice answer. I would like to provide another one.

Consider a charge $q$ be enclosed by any surface (not necessarily a sphere). Something like this -

Now, you write the flux coming out of this weird surface - $$ \phi_E = \displaystyle \oint_S \mathbf{E} \cdot \mathrm{d\mathbf{S}} $$ We know that - $$ \mathbf{E} = E \vec{r} = \dfrac{1}{4\pi \epsilon_0} \dfrac{q}{r^2}\vec{r} $$ So, here in this weird surface. there is no fixed radius, is there? And the surface area that is considered here is not a continuous one. So, I would get - $$ \phi_E = \dfrac{q}{4\pi \epsilon_0} \displaystyle \oint_S \dfrac{\mathrm{d\mathbf{S}}}{r^2} \tag{1} $$ Recall that the term $\dfrac{\mathrm{d\mathbf{S}}}{r^2}$ is the very definition for the steradian - which is equal to $\dfrac{1}{4\pi}$ of a complete sphere. This holds good for any surface. Simply put, this is the 3D analogue of the $2\pi$ rotation in a circle. Here, we have its differential element, i.e, $d\Omega = \dfrac{\mathrm{d\mathbf{S}}}{r^2}$ Integrating it entirely, we have $$ \displaystyle \oint_S \dfrac{\mathrm{d\mathbf{S}}}{r^2} = \displaystyle \oint_S d\Omega = 4\pi$$ Plugging this back into (1), we have - $$ \phi_E = \dfrac{q}{\epsilon_0} $$ Which implies - $$ \displaystyle \oint_S \mathbf{E} \cdot \mathrm{d\mathbf{S}} = \dfrac{q}{\epsilon_0} $$ Ok, since we're done with deriving the integral form of Gauss's law (Which holds true for any closed surface), the following differential form can be obtained by applying the divergence theorem - $$ \nabla \cdot \mathbf{E} = \dfrac{\rho}{\epsilon_0} $$