Then it goes on explaining how Gauss law would fail because for a very large surface, E field would be vanish with flux through it and though we can calculate div for this field it won't depend on source density.

But I don't get what makes this particular function so evil that it would break physics, I can see that fallout rate is exponentially large for this function than would for a "pure inverse relationship".

Basically I want to know what makes this function so different than the ##1/r^2## relationship ?

"Break physics" is too vague. I think the main point is that there are consequences beyond the obvious. Maxwell's equations and all of electromagnetism would be much different.

Physics would still work, but it would work much differently from what we are used to.

Yes, but things that fall off faster than ##1/r^2## turn out to be much different from things that fall off slower than ##1/r^2##.

##1/r^2## is unique and special for lots of reasons, among them are Gauss's law and Maxwell's equations. But there are also some symmetries to ##1/r^2## that give rise to special and unique things in quantum mechanics: for example the degenerate energy levels of the hydrogen atom. Atomic physics and chemistry will all change significantly if there are significant deviations from the existing Coulomb's law.

The Yukawa potential is the static potential due to a field with mass. In this case instead of the Poisson Equation you get
$$(\Delta+m^2)\phi=-\rho.$$
Now to get the Green function, i.e., the solution for a poin charge ##\rho=Q \delta^{(3)}(\vec{x})##, use a Fourier transform,
$$\phi(\vec{x})=\int_{\mathbb{R}^3} \mathrm{d}^3 \vec{k} \frac{1}{(2 \pi)^3} \tilde{\phi}(\vec{k}) \exp(\mathrm{i} \vec{x} \cdot \vec{k}).$$
Inserting this ansatz into the equation you get
$$(k^2+m^2) \tilde{\phi}=1 \; \Rightarrow \; \tilde{\phi}(\vec{k})=\frac{1}{k^2+m^2}.$$
For the Fourier integral you introduce spherical coordinates with the polar axis along ##\vec{x}##. After some algebra you are left with
$$\phi(\vec{x})=\frac{1}{4 \pi^2 r} \int_{-\infty}^{\infty} \mathrm{d} k \frac{k}{k^2+m^2} \sin(k r).$$
This integral can be solved with the usual trick making ##k## complex and close the contour by a large semicircle in the upper (lower) plane for the ##\exp(\mathrm{i} k r)## (##\exp(-\mathrm{i} k r)##) part of the sine. You pick up the poles at ##\mathrm{i} m## (##-\mathrm{i}m##), respectively, leading finally to
$$\phi(\vec{x})=\frac{1}{4 \pi r^2} \exp(-m r).$$
Of course, for ##m=0## you get back the Green's function for ##-\Delta##, i.e., the Coulomb potential for a unit charge.