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
There have been lots of experimental attempts to test the validity of Coulomb's $r^{-2}$ law. Many of these are reviewed by Tu & Luo (2004), and is where I am getting the numbers quoted below. Somewhat equivalently, experiments have looked at trying to set an upper limit to the photon mass, which is testing the hypothesis that rather than a $r^{-1}$ relation, that the Coulomb potential falls in a similar way to the Yukawa potential, as $r^{-1} \exp(-m_\gamma c r/\hbar)$.
The laboratory tests largely involve measuring the potentials on concentric charged spheres and are relatively small scale. These show that if Coulomb's law scaled as $1/r^{2+q}$, then the current limits are $|q|< 10^{-16}$. On the (laboratory) scales probed by the experiments, this corresponds to an upper limit to the photon mass of $m_\gamma < 10^{-50}$ kg (Crandall 1983; Fulcher 1986).
The size of laboratory equipment limits the constraints one can put on the mass of the photon and the scale-length of any Yukawa-like potential. However, on large scales, a non-zero photon rest mass would lead to a number of observational effects. Not only is the potential changed, but there is a predicted frequency-dependent velocity and the possibility of longitudinally polarised photons. The most stringent limit appears to come from considering the stability of magnetised gas in galaxies, where the claim is that the photon mass must be less than $10^{-62}$ kg, which is equivalent to a Yukawa-like scale length of 1000 pc! (Chibisov 1976). It is not clear how seriously this claim is taken, but Tu & Luo (2004) list several other cosmological and laboratory studies that have placed limits on any scalelength of $>10^{10}$ m. At a distance of 1000 km, these deviations would amount to a difference in force of $\exp(-1000)$.
So from the point of view of your question, there is experimental evidence that the deviations from the Coulomb law are utterly negligible at scales of 1000 km.