I am not quite sure which literature are you looking for, but it should be written in standard textbook. Anyway, you are right that it is not that useful except a formulation and exercise.
Usually, it is written in the following form. Given a charge distribution $\rho_{1}(\vec{\mathbf{r}_{1}})$ for object with volume $V_1$ and $\rho_{2}(\vec{\mathbf{r}_{2}})$ for object with volume $V_2$, we have the electric field:
$$\vec{\mathbf{E}}(\vec{\mathbf{r}_{2}})=\frac{1}{4\pi\epsilon_{0}}\int_{V_{1}}\frac{\rho_{1}(\vec{\mathbf{r}_{1}})}{|\vec{\mathbf{r}_{2}}-\vec{\mathbf{r}_{1}}|^{3}}(\vec{\mathbf{r}_{2}}-\vec{\mathbf{r}_{1}})d\vec{\mathbf{r}_{1}} \tag{1}$$
and force
$$\vec{\mathbf{F}}=\frac{1}{4\pi\epsilon_{0}}\int_{V_{2}}\int_{V_{1}}\frac{\rho_{1}(\vec{\mathbf{r}_{1}})\rho_{2}(\vec{\mathbf{r}_{2}})}{|\vec{\mathbf{r}_{2}}-\vec{\mathbf{r}_{1}}|^{3}}(\vec{\mathbf{r}_{2}}-\vec{\mathbf{r}_{1}})d\vec{\mathbf{r}_{1}}d\vec{\mathbf{r}_{2}} \tag{2}$$
Because of practicality, it is not used that often:
- It is hard to apply on non-fixed charge distribution. For example, for conductor, the charge itself would change as two object approaching. This internal dynamics is not captured in Eq (2).
- When the objects are far away, it is just like a point charge. Why do you want to do the cumbersome integration?
- There are better tools to handle this situation: Multipole expansion. As two objects become closer and closer, you include the monopole first (i.e. a point charge), then dipolar, then quadrupole... This expansion is very systematical and have good physical meaning.
- People care more about electric field rather than force, as it is more fundamental, so the Eq (1) is more emphasized than Eq (2).
I believe the problem with your proof is that packing an area with spheres only takes up about 75% of a volume, and this problem doesn't disappear when you shrink the spheres down to infinitesimal sizes. No matter how small you shrink the spheres, those same packing properties should hold.
Your proof relies on saying that the integral over a large surface is equivalent to adding up tiny infinitesimal spherical surfaces. But packing with spheres won't give 100% efficiency. So the volume of the small spheres won't add up to the volume inside of the larger surface, and the two integrals $\int V d\tau$ won't be equivalent because they span different volumes.
There's the same kind of problem with the surface integral. If you try to add two shapes together with a flush side, the electric field will be equal while the normals will be in opposite directions, so the flush surface cancels out and you're left with a surface integral of just the boundary of the two shapes. With spherical packing, the surfaces will not be flush, and so the surface integral of the larger shape will be different than for the smaller shape.
The book proves you can make a dent in the sphere and the integral should remain the same. Therefore, you can draw any shape you want around a point charge and the integral should remain the same. This is equivalent to saying you can move a point charge anywhere you want inside an arbitrary surface. Now you can use superposition. If the law holds for any point charge at any position inside a surface, then break up your arbitrary charge distribution into infinitesimal point charges $dq$. The law holds for all $dq$, so it should hold for $\int_V dq d\tau$.
Proving a 'dent' in a sphere doesn't change the integral is equivalent to proving that for some arbitrary volume element, the integral is zero. This is because making a dent is equivalent to subtracting a small volume element. This is actually not too difficult. We can find $V$ just by inspection using $V=-\nabla E$.
$$ V= \frac{q}{4\pi \epsilon_0 k}e^{\frac{-k}{\lambda}}$$
Convert the surface integral of the electric field into a volume integral with the divergence theorem. The field only depends on $r$ so using the spherical form of $\nabla$ we get
$$ \nabla \cdot \vec{E }= \frac{1}{k^2}\frac{\partial}{\partial k} (k^2\vec{E})= \frac{1}{k^2}\frac{q}{4\pi \epsilon_0}\left( \frac{1}{\lambda} + (1 +\frac{k}{\lambda})\frac{-1}{\lambda}\right)e^{\frac{-k}{\lambda}}$$
$$ = \frac{q}{k^24\pi \epsilon_0} \frac{-k}{\lambda^2}e^{\frac{-k}{\lambda}}$$.
Now simply rewrite the new Gauss's law.
$$ \int_V \nabla \cdot \vec{E} d\tau + \frac{1}{\lambda^2}\int_V V d\tau= \frac{q}{4\pi \epsilon_0 \lambda^2}\int_V \left( \frac{-1}{k}e^{\frac{-k}{\lambda}} + \frac{1}{k}e^{\frac{-k}{\lambda}}\right)d\tau.$$
At $k=0$ the equation blows up on account of the point charge. For other finite $k$ values, this equation is just zero. This means that for any volume element that does not include the point charge, the two integrals add to zero. This means you can add or subtract any crazy surfaces you want from your original sphere as long as they don't include the origin.
Best Answer
He used what is called a torsion balance. His experimental method is outlined very nicely in this video.
After Coulomb published the result of his work, however, it was debated as to whether his experiment really did provide enough evidence to support his claim that the force between two point charges really did follow the equation we now call Coulomb's Law. This was because subsequent attempts to recreate his experiment proved that it was very difficult to obtain consistent results using Coulomb's method. If you are interested, here (DOI) is a very detailed report which outlines one of these recreations. This author comes to the conclusion that Coulomb's experiment was very very sensitive to any slight change: "Thus we see how just a single, slight, and easily invisible difference between Coulomb’s apt prescriptions and a replication may lead to entirely different results" (561).
Of course, it was later shown more solidly that Coulomb's law does hold up to experimental evidence, but it is still debated as to whether or not Coulomb jumped to conclusions, so to speak, and simply got lucky.