Coulomb’s Law
Consider the field of a single positive point charge q and place this charge at the center of an imaginary spherical surface S, with radius r. From Coulomb’s law, we know that:
$$\vec{E} = \frac{1}{4 \pi \epsilon_{0}} \frac{q}{r^{2}} \, \hat{r}$$
The magnitude of electric field is then:
$$E = \frac{1}{4 \pi \epsilon_{0}} \frac{q}{r^{2}}$$
We will now attempt to get the essence of Gauss’s law from Coulomb’s law. Rearranging the above equation:
$$E \times 4 \pi r^{2} = \frac{q}{\epsilon{0}}$$
We will focus on the left hand side of the equation first. Note that $4 \pi r^{2}$ is the surface area of a sphere.
$$E \times 4 \pi r^{2} = E \times \oint\limits_{S_{r}} dA$$
Since E is a constant everywhere on the surface of $S_{r}$, we can bring E into the integral.
$$\begin{aligned} E \times 4 \pi r^{2} &= \oint\limits_{S_{r}} E \, dA \\ &= \oint\limits_{S_{r}} \vec{E} . d \vec{A} \end{aligned}$$
Recall that $d \Phi_{E} = \vec{E} . d \vec{A}$. Hence, the above equation is just $\Phi_{E}$. Combining with the right hand side of the original equation from Coulomb’s Law, you obtain $\Phi_{E} = \frac{q}{\epsilon_{0}}$.
This says that the total electric flux through the spherical surface is given by $\frac{q}{\epsilon_{0}}$, where q is the charge enclosed by the spherical surface.
Notice that the total electric flux is independent of the radius r of the sphere. We can see this by calculating the electric flux for two spheres of radius R and 2 R respectively.
For sphere of radius R:
$$ E = \frac{1}{4 \pi \epsilon_{0}} \frac{q}{R^{2}}$$
$$\begin{aligned} \Phi_{E} &= \oint\limits_{S_{R}} \vec{E}.d\vec{A} \\ &= \oint\limits_{S_{R}} E \, dA \\ &= E \oint\limits_{S_{R}} \, dA \\ &= E \times 4 \pi R^{2} \\ &= \frac{q}{\epsilon_{0}} \end{aligned}$$
For sphere of radius 2R:
$$E = \frac{1}{4 \pi \epsilon_{0}} \frac{q}{(2R)^{2}}$$
$$\begin{aligned} \Phi_{E} &= \oint\limits_{S_{2R}} \vec{E}.d\vec{A} \\ &= E \times 4 \pi (2R)^{2} \\ &= \frac{q}{\epsilon_{0}} \end{aligned}$$
You can continue on for all the different radius, and this result will still hold. This is because as the radius increases, E decreases by a factor of $\frac{1}{r^{2}}$. In other words, this result will only hold if the $E \propto \frac{1}{r^{2}}$ exactly.
Gauss’s law will hold for a surface of any shape or size, provided that it is a closed surface enclosing the charge q.
Note: We have “shown” that Gauss’s law is compatible with Coulomb’s law for spherical surfaces. What about non-spherical surfaces? I am afraid that you will have to take my word for it. You will need the ideas of “divergence” to have a more proper proof.
