I really appreciate the physical explanations made in other answers, but I want to add that Fourier transform of the Coulomb potential makes mathematical sense, too. This answer is meant to clarify the sense that the standard calculation is mathematically valid.
Firstly, and maybe more importantly, I want to emphasize that
The Fourier transform of $f$ is not simply just $\int{f(x)e^{-ikx}dx}$.
For an $L^1$ function (a function that is norm integrable), this is always the case but the Coulomb potential is definitely not in $L^1$. So the Fourier transform of it, if it ever exists, is not expected to be the integral above.
So here comes the second question: Can Fourier transformation be defined on functions other than $L^1$?
The answer is "yes", and there are many Fourier transformations. Here are two examples.
- Fourier transformation on $L^2$ functions (i.e., square integrable functions).
It turns out that the Fourier transform behaves more nicely on $L^2$ than on $L^1$, thanks to Plancherel's theorem. However, as we mentioned above, if an $L^2$ function is not in $L^1$, then the above integral may not exist and the Fourier transform is not given by that integral, either. (However, it has a simple characterization theorem, saying that in this case the Fourier transform is given by the principle-value integration of the above integral.)
- Fourier transform of distributions (generalized functions)
It is in this sense that the Fourier transform of Coulomb potential holds. The Coulomb potential, although not an $L^1$ or $L^2$ function, is a distribution. So we need to use the definition of the Fourier transform to distributions in this case. Indeed, one can check the definition and directly calculate the Fourier transform of it. However, the physicists' calculation illustrates another point.
Fourier transformation on distributions (however it is defined) is continuous (under a certain topology on the distribution space, but let's not be too specific about it).
Remember that if $f$ is continuous, then $x_\epsilon\rightarrow x$ implies that $f(x_\epsilon)\rightarrow f(x)$.
Now $\frac{1}{r}e^{-\mu r}\rightarrow\frac{1}{r}$ when $\mu\rightarrow 0$ (again, under the "certain topology" mentioned above), and therefore continuity implies
$$\operatorname{Fourier}\left\{\frac{1}{r}e^{-\mu r}\right\}\rightarrow \operatorname{Fourier}\left\{\frac{1}{r}\right\}.$$
However, $\frac{1}{r}e^{-\mu r}$ is in $L^1$ and therefore its Fourier transformation can be computed using the integral $\int{f(x)e^{-ikx}dx}$.
Therefore those physicists' computations make perfect mathematical sense, but it's on Fourier transform of distributions, which is much more general than that on $L^1$ functions.
I hope this answer can build people's confidence that the Coulomb potential Fourier transformation problem is not only physically reasonable but also mathematically justifiable.