For proving that the Coulomb E satisfies the Gauss law in differential form, simply take the divergence of Coulombian E. You will need the Dirac delta and some of its properties.
For deriving the divergence of E field, do just the same as above...
You can also derive the Gauss law in integral form from the differential form by applying the divergence theorem (= Gauss theorem --- I guess you know, why :)).
The integral form of the Gauss law from the Coulomb E can be directly derived by writing the flux integral for the Coulomb E, using the linear independency of the two integrals (primed vs. non-primed), and then recognizing a common vector integral expression.
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
The infinitely long wire has an infinite charge $$Q~=~\lambda \int_{-\infty}^{\infty} \! dz ~=~ \infty,\tag{1}$$ and EM has an infinite range, so one shouldn't be surprised to learn that the result $$\phi(r)~=~ \frac{\lambda}{4 \pi \epsilon_0} \int_{-\infty}^{\infty} \frac{dz}{\sqrt{z^2+r^2}} ~=~ \frac{\lambda}{4 \pi \epsilon_0} \left[ {\rm arsinh} \left(\frac{z}{r}\right)\right]_{z=-\infty}^{z=\infty} ~=~\infty \tag{2}$$ is infinite. (From a mathematical point of view, the integrand fails to be integrable wrt. the $z$ variable.) A similar situation happens often in Newtonian gravity if the total mass is infinite, see e.g. this question.
However, as Pygmalion mentions in his answer, the electric field $\vec{E}$ is well-defined for $r\neq 0$, and the corresponding integrand is integrable wrt. the $z$ variable. E.g. the radial component (in cylindrical coordinates) reads $$E_r(r)~=~ \frac{\lambda r}{4 \pi \epsilon_0} \int_{-\infty}^{\infty} \frac{dz}{(z^2+r^2)^{3/2}} ~=~ \frac{\lambda r}{4 \pi \epsilon_0} \left[ \frac{z}{r^2\sqrt{z^2+r^2}}\right]_{z=-\infty}^{z=\infty} ~=~\frac{\lambda}{2\pi\epsilon_0 r}\tag{3} $$ for $r\neq 0$.
Alternatively, apply Gauss' law $$d\Phi_{E} ~=~\frac{dQ}{\epsilon_0},\tag{4}$$ using an infinitesimally thin disk perpendicular to the wire. The disk has radius $r$ and thickness $dz$. The total electric flux $d\Phi_{E}$ out of the disk is $$ E_r \cdot 2\pi r dz ~=~ \frac{\lambda dz}{\epsilon_0},\tag{5}$$ which leads to the same electric field $E_r$.
This electric field $\vec{E}=-\vec{\nabla}\phi$ is consistent with a potential of the form $$\phi(r) ~=~-\frac{\lambda}{2\pi\epsilon_0}\ln r \qquad \text{for}\qquad r\neq 0.\tag{6}$$