I wondered about this question for a while myself, but I think by now I've become convinced why there really are so many smooth functions. Maybe the ideas that convinced me can convince you. I think you answered the question when you asked it: the mollification process is basically responsible. If I recall correctly, the business about smooth homotopies is also most easily proven via mollification.
If you start with a non-negative smooth bump function $\eta$ in the unit ball of ${\mathbb R}^n$ normalized so that $\int \eta(x) dx = 1$, then you can regard the measure $\eta(x) dx$ as a smooth probability measure. Likewise the measure $\eta_\epsilon(x) dx = \eta(\frac{x}{\epsilon}) \frac{dx}{\epsilon^n}$ is another smooth, probability measure supported in the ball of radius $\epsilon$ (it's the pushforward of the first one by multiplication by $\epsilon$).
So with this interpretation, the mollification $f_\epsilon(x) = \int f(y) \eta_\epsilon(x - y) dy$ is what happens when you randomly translate the function $f$ and replace its value at each point by the expected value after these random translations. Now, it should seem intuitive that even when the original function is singular, the resulting function looks a lot like the original one, but is a whole lot smoother because the singularities get spread out and hence diminished. Imagine, for example, what results when you do this to the characteristic function of an open set -- you get a smooth cutoff that looks like the rough one. Or if you want a smooth partition of unity, start with a rough one (characteristic functions of sets) and just mollify it to get a smooth one.
So, if you believe in smooth bump functions you should believe there are a lot of smooth functions. By the way, in case you don't know this, there is more than one way to produce a smooth bump function. One way (say we're on the line) is to repeatedly convolve characteristic functions of intervals in such a way that the support remains bounded -- the regularity increases by one every time you do it.
EDIT: I noticed that the question also asked about what makes smooth functions so useful, which is something I didn't address at all.
One reason they are so useful is that they simply do not have the defects that non-smooth functions have, and so they are less likely to introduce issues which are irrelevant to your problem. For example, if you want an asymptotic count of the integer lattice points in a large ball, you want to use Poisson summation so that the main term is simply the contribution from the $0$ frequency and everything else can be treated as an error. Unfortunately, the Fourier transform of the characteristic function of a ball fails to decay well enough for this idea to work (although it decays better than you would expect thanks to the curvature of the sphere, a fact which ultimately improves the error term). The problem here is related to the uncertainty principle -- Fourier analysis cannot give a count of lattice points that is precisely sensitive to points along the boundary; the characteristic function of an open ball or its closure define the same $L^2$ function. Thus, you have to mollify the ball for this strategy to work. Essentially the same issue arises in proving the prime number theorem (and commonly arises in analytic number theory in general) -- many proofs have irrelevant technicalities, all related to failing to smooth out the count.
Another instance: proving that the fundamental group of the 2-sphere is trivial. It's easy to show that a smooth (or Lipschitz) curve will miss some point on the $2$ sphere (since it can't even increase Hausdorff dimension), and therefore such a curve is homotopic to a constant. But continuous curves could cover the whole sphere, giving a problem related more to oscillations and a lack of regularity and you might say less to topology. A similar problem arises when considering the degree of a continuous map -- they can hit their values infinitely many times thanks to uncontrolled oscillations so it's harder to interpret the degree as a count and usually requires an approximation. On the other hand, the distinctions between topological and smooth manifolds (e.g. exotic spheres) makes a really interesting topic, so one can't always dismiss lack of smoothness as just being a nuisance.
One more unifying reason smooth functions are so useful is the fact that singular functions can be approximated by smooth functions, so concepts which already have explicit meaning for smooth functions (e.g. degree, distribution theoretic derivatives, Fourier transform) have a natural extension to singular functions exactly when they are continuous with respect to that kind of approximation. For example, degree is continuous with respect to uniform approximation so continuous functions have a degree (more interesting is that degree extends to maps with approximations in BMO -- see surveys of Brezis). When you work with smooth functions you rarely worry that something is defined, and then you just have to check the estimates/continuity (e.g. Fourier Transform maps $L^2 \to L^2$) to make sure defining by approximation is OK.
Finally, you can only use smooth functions when you are on a manifold, so many of the axioms that people write down in point set topology to make sure they are not looking at some atrocious, pathological space, have already been built (entirely?) into the machinery of smooth functions and partitions of unity. Smooth functions can tell the difference between two closed sets, therefore the space is normal, etc. Thus, using smooth functions in arguments helps you avoid unwanted pathologies of the space as well as the maps.