Gauss' law for two point charges and a factor of $4\pi$ When we model electric field as something that represents the "flow" of something that is conserved, we can prove that the flux due to a single point charge, through the following surface, 
Flux through a radial surface
is effectively zero. And therefore prove (almost) that flux through all kinds of  surfaces that do not include a charge is zero. But when we do have a charge, we perform a clever trick, we enclose the charge with another surface, in the following manner, 
Flux due to point charge
and therefore prove that the flux through the required surface is equal to the flux through small surface enclosing the charge.
Therefore, 
$$ \int_S \mathbf E . \mathbf {ds} = \dfrac{q}{4\pi \epsilon r^2} (4\pi r^2) \\ \int_S \mathbf E . \mathbf {ds} = \dfrac{q}{ \epsilon} 
 $$
This "trick" was possible, only because,  we had a point charge and the field was radially outward. So my question is, how are we able to generalise that the flux through any surface is equal to this fraction $\dfrac{q}{\epsilon}$? What if we had two point charges instead of one, then we wouldn't get a convenient $4\pi r^2$ term.  This again leads to a question; why in the first place, is there a $4\pi$  in Coulomb's law? 
P.S: I'm relatively new to the subject, so I'm missing something obvious.And I couldn't choose a suitable title for the question:|
 A: Why does Gauss' law (in integral form) hold for any closed surface? One way of looking at it is that it just does: it's one of Maxwell's equations, and there is no "proving" it. In other words, there is experimental evidence for it. Many people find it more natural to express it in the differential form ($\vec{\nabla}·\vec{E} = \rho/\epsilon_0$), which is equivalent to the integral form (through the divergence theorem) but doesn't refer to a specific surface. Gauss' law can be motivated by Coulomb's law, or even derived from it in electrostatics, but Gauss' law is more fundamental in the sense that it holds in electrodynamics while Coulomb's law does not. 
Here is a somewhat lazy derivation of Gauss' law from Coulomb's law using the divergence theorem that also illustrates how $4\pi$ comes about. In general it falls out of the math when dealing with situations of spherical symmetry.
From Coulomb's law,
$$\vec{E}(\vec{r})=\frac{1}{4\pi\epsilon_0}\iiint\limits_\infty dV' \rho(\vec{r}') \frac{\vec{r} - \vec{r}'}{\lVert \vec{r} - \vec{r}' \rVert^3}.$$
Consider a region $V$ bounded by the closed surface $\partial V$. Then,
$$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iint\limits_{\partial V}{ \iiint\limits_\infty dV' \rho(\vec{r}') \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3} }.$$
Changing the order of the integrals (which I won't attempt to justify),
$$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iiint\limits_{\infty}{ dV' \rho(\vec{r}') \iint\limits_{\partial V} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3} } }$$
Now consider the integral
$$ I(\vec{r}') = \iint\limits_{\partial V} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}} .$$
Suppose that $\vec{r}'$ is not enclosed by $\partial V$, so that the denominator is always non-zero and the integrand is always well defined. By the divergence theorem,
$$ I(\vec{r}') = \iiint\limits_V {dV \vec{\nabla} · \frac{\vec{r} - \vec{r}'}{\lVert \vec{r} - \vec{r}' \rVert^3}} . $$
It's a straightforward calculation to show that the divergence in the integrand vanishes, therefore $I = 0$.
If on the other hand $\partial V$ encloses $\vec{r}'$, let $S_\delta$ be a sphere of radius $\delta$ centered about $\vec{r}'$ and contained in $\partial V$. Let $V - S_\delta$ be the subset of $V$ outside $S_\delta$, $\partial(V - S_\delta)$ the closed boundary of this region and $\partial S_\delta$ the spherical boundary of $S_\delta$. We can write
$$ I(\vec{r}') = \iint\limits_{\partial(V - S_\delta)} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}} + \iint\limits_{\partial S_\delta} { \frac{(\vec{r} - \vec{r}')·d\vec{S}}{\lVert \vec{r} - \vec{r}' \rVert^3}}. $$
The first integral is zero for the same reason as before: the divergence of the integrand vanishes everywhere. The second integral can now be evaluated easily using spherical symmetry to give
$$ I(\vec{r}') = \frac{\delta\times4\pi\delta^2}{\delta^3} = 4\pi.$$ This is where the factor of $4\pi$ comes from. Thus,
$$ \iint\limits_{\partial V}{\vec{E}(\vec{r})} · d\vec{S} =\frac{1}{4\pi\epsilon_0} \iiint\limits_{\infty}{ dV' \rho(\vec{r}') I(\vec{r}) } = \frac{1}{\epsilon_0} \iiint\limits_{V}{ dV' \rho(\vec{r}') } = \frac{Q}{\epsilon_0} $$
where the second equality uses the fact the $I(\vec{r}')$ is non-zero if and only if $V$ contains $\vec{r}'$.
A: The key point of the argument is that the net flux does not depend on the shape of the surface enclosing the charges.   This is a physical argument, not a mathematical one.  
Thus, if you have a net flux through a small spherical surface, and this small surface is inside a box, then the net flux through the box is the same as the flux through the spherical surface (barring sources or sinks elsewhere in the box), simply because every field line that escapes the sphere must also escape the box.
Now, if you were to compute the flux through the sides of a box rather than a sphere, the flux would NOT be uniform through every surface element of the sides.  Thus, the integral 
$$
\int \vec E\cdot d\vec S
$$
is not easy to evaluate because, as $\vec E$ is not constant over the surface, and because the angle between $\vec E$ and $d\vec S$ is not constant over the surface; indeed this integral does not equal $\vert \vec E\vert S$ so there’s no $4\pi$ factor in the first place (which $\vert \vec E\vert$ would you use anyways since it is not constant over the surface?).  You would still know the net flux - it *must * be $q/\epsilon_0$ - but you cannot say anything more.  
This is to be constrasted with the spherical surface.  By symmetry, the field has constant magnitude on the surface and also must point everywhere radially.  Thus, for this surface and this surface only, 
\begin{align}
\int \vec E\cdot d\vec S = \int \vert \vec E\vert dS = \vert \vec E\vert \int dS
=4\pi r^2 \vert \vec E\vert
\end{align}
where, in going from the second to the third integral, the observation that 
$\vert \vec E\vert$ is constant on the surface has been used to pull it out of the integral.  
In a nutshell, the computation of the flux by integration can only be done in highly symmetrical cases where $\vert \vec E\vert$ is constant on the surface, but the net flux does not depend on this symmetry because physically it’s a function only of sources and sinks inside the Gaussian surface.  Thus, to compute the net flux one choose one (or more) convenient surface to evaluate the integral, then invokes the physics to extend to any surface enclosing the charges.
A: Put a point charge inside a surface of any shape (or size). Consider  a thin cone which extends out from the charge. Where the cone intersects the surface, the surface area is represented by a vector pointing outward.  The associated flux is defined as the dot product of the E field and the area vector.  We want the component of the field which is crossing the area, but that can also be interpreted as using the component of the area vector which is parallel to the field.  This component of the area increases with the square of the distance as you move out along the cone, and the field decreases with the square of the distance.  The flux out through the cone is a constant.  Combine all the cones coming out from the charge and you get a constant total flux, which is independent of the size or shape of the surrounding surface.  Choose a surface of constant radius, and the flux becomes E(4pi R^2).  Each additional charge contributes its own flux.
