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.
This is all about the structure of the equations. The electrostatic field $\mathbf{E}$ satisfies the two equations
$$\nabla\cdot \mathbf{E}=\dfrac{\rho}{\epsilon_0},\quad \nabla\times \mathbf{E}=0,$$
together with boundary conditions on crossing charged surfaces, for example.
What the author means is the following. We know that $\mathbf{D}$ satisfies the equations
$$\nabla\cdot \mathbf{D}=\rho_{\text{free}},\quad\nabla\times\mathbf{D}=\nabla\times \mathbf{P}.$$
Now when $\nabla\times \mathbf{P}=0$, for example, when $\mathbf{P}$ is constant, we have that
$$\nabla \cdot \mathbf{D}=\rho_{\text{free}},\quad \nabla \times \mathbf{D}=0.$$
As a mathematical problem of solving a differential equation this is the same as the electrostatic problem I said above, with the identification that the charge density should be $\rho = \epsilon_0\rho_{\text{free}}$.
To see this, divide by $\epsilon_0$ the first equation and notice that
$$\nabla\cdot \left(\dfrac{1}{\epsilon_0}\mathbf{D}\right)=\dfrac{\rho_\text{free}}{\epsilon_0},\quad \nabla \times \mathbf{D}=0.$$
Thus since these are the same equations for the electrostatic field in vacuum generated by a charge density $\rho_{\text{free}}$, which we shall call $\mathbf{E}_{\text{vac}}$ you have that
$$\dfrac{1}{\epsilon_0}\mathbf{D}=\mathbf{E}_{\text{vac}}$$
which gives your equation $\mathbf{D}=\epsilon_0 \mathbf{E}_{\text{vac}}$.
Best Answer
Coulomb did not know the absolute value of the charge but what he was able to do was to reduce the charge on one of his spheres by a known ratio.
He charged a metal sphere and used it in his experiment.
He then removed that metal sphere and touched it with an identical uncharged metal sphere.
He assumed that the final charge of the initially charged sphere was half of what had been initially because the initially uncharged sphere has removed half the charge from the initially charged sphere.
He could then use that initially charged sphere with half its initial charge to take a second set of readings.
You should note that defining the coulomb as the unit of charge is a rather recent.
The electrostatic unit (esu) of charge was defined using Coulomb's law, $F = \dfrac {q_1\,q_2}{r^2}$, where $F$ is the force of attraction/repulsion in dynes (the force required at accelerate $1$ gramme at $1\, \rm cm \,s^{-2}$), $r$ the separation in centimetres and the charge was then in esu or statCoulomb or franklin.
There have been many who have tried to reproduce Coulomb's original experiments to try and evaluate the sort of accuracy that Coulomb might have been able to achieve and even as to whether Coulomb actually got his "results" experimentally.
The paper The Material Intricacies of Coulomb’s 1785 Electric Torsion Balance Experiment and the links therein may be of interest?