Why are there so many smooth functions? I do understand that my question might seem a little bit ignorant, but I thought about it a lot and still can't wrap my head around it. 
Analycity imposes very strong conditions on a map, from elementary ones like "locally zero implies globally zero", to a little bit more deep like the Hurwitz formula (in the complex case). Neither of above are true if we just assume smoothness. 
On one hand, it is quite easy to prove that smooth functions are dense in any "reasonable" function space (I guess it depends on what one considers reasonable, though…) - just convolve with smooth approximations of identity. Also, (although I do take it on faith), any map of two manifolds is homotopic to a smooth one and two homotopic smooth maps are actually smooth-homotopic. 
Because of above facts, it seems to me that smooth functions are "abundant" and are actually very close to topology, ie. mere continuity. 
On the other hand, when I recall basic calculus course, it always seemed like being differentiable even once is a "miracle", and being differentiable infinitely many times is a very, very strong condition, even more so in several variables. 
Why are objects so constrained, ie. smooth functions, so useful and also, so malleable?
 A: Smoothness only constrains the infinitesimal behaviour of a function $f$ at a given point $x_0$ - the limiting behaviour of $f(x)$ as $x \to x_0$.  Analyticity constrains the local behaviour - the value of $f$ on a non-infinitesimal ball $B(x_0,\varepsilon)$ - a far stronger constraint.  
Actually, one should make a distinction between qualitative smoothness - i.e. infinite differentiability, with no bounds on the derivatives - and quantitative smoothness - things like bounds on $C^k$ norms (or variants of these norms, such as Holder norms or Sobolev norms).  The latter does control local behaviour and not just infinitesimal behaviour, thanks to tools such as Taylor's theorem with remainder, and is a far stronger constraint; a highly oscillatory function need not be close in reasonable norms to smooth functions with good bounds on $C^k$ norms.
The complex variable case is a bit different from the real variable case, as the mere fact of differentiability now imposes an elliptic constraint, namely the Cauchy-Riemann equations: see Why do functions in complex analysis behave so well? (as opposed to functions in real analysis)
A: On the other hand, in the Banach space $C[0,1]$ with the supremum norm, the set of functions which do not have (finite) derivative at any of the points of $[0,1]$ is residual. This was first proved by Stefan Banach and Stefan Mazurkiewicz in 1930s. A proof following J. C. Oxtoby can be found here (if you read Spanish): 
http://www.ciens.ula.ve/matematica/publicaciones/libros/por_profesor/wilman_brito/Teorema_de_Baire_Aplicaciones.pdf
For an even more refined result (a combination of several ones), look up Banach-Mazurkiewicz-Jarnik theorem (Theorem 7.2.1) in the following book:
Marek Jarnicki, Peter Pflug: Continuous nowhere differentiable functions. Monsters of analysis. Springer Monographs in Mathematics, 2015.
A: One way to think about this is through solutions of (funcional) equations. If you are studying solutions to polynomial equations, then you are in the realm of algebraic geometry, and similarly there is the study of solutions to  equations of analytic functions. The depth of these fields properly justifies the claim that we have made "strong restrictions" on the functions we allow. 
However, in this context, smoothness imposes no such restrictions. In fact many theorems on the continuous setting extend to the smooth setting (with a bit of machinery of course, partitions of unity or bump functions), for example one has a smooth Urysohn lemma which implies that studying zero sets of smooth functions is the same as studying closed sets of our space.
From the most basic examples one also notices that for "nice enough" functional equations imposing smoothness on the solutions is the same as imposing conitnuity. The first time one probably encounters this is in Cauchy's functional equations and its analogues, which say that the solutions are either smooth (linear, polynomial, exponential, logarithmic etc.) or they are wild (dense in the plane). I don't know of a way to make this functional equation business rigorous in general, but in case one restricts to group structures then this phenomenon is basically Hilbert's fifth problem, a topological group which is locally Euclidean is a Lie group. Here is a nice exposition by Terence Tao.
A: Here is a somewhat different answer:
Smooth functions always admit asymptotic expansions in each point, convergent or not. One difference between analytic functions and smooth functions is, of course, that analytic functions are uniquely determined by their asymptotic expansion at some point (which, of course, is just their Taylor series), while smooth functions are not. But even if smooth functions were determined by their asymptotic expansion at a given point, you still see what goes wrong:
There are Theorems (as far as I recall by Witt for the complex case and by Borel for the real-analytic case, though both only in one dimension) that answer the question what constraints a sequence has to fulfill to be the coefficient series of the asymptotic expansion of a smooth function. The answer: There is none.
So, every real or complex sequence appears as coefficient sequence of the asymptotic expansion of a smooth function at a given point! 
This means, the space of possible coefficients is bigger than $L^{\infty}(\mathbb{N},\mu)$ (for any measure equivalent to the point measure) and that space isn't separable. And it gets even worse because a given coefficient series doesn't locally determine a smooth function, much less globally. Coefficient series of an analytic function however lie in the space $L^1(\mathbb{N},\mu)$, where $\mu$ is the point measure multiplied with the function $n \mapsto 1/n!$, which is separable.
A: Given a paracompact smooth manifold, you have smooth partitions of unity (nLab), but on a real analytic manifold (e.g. a complex manifold viewed as a real manifold) one doesn't have analytic partitions of unity (much less holomorphic, if you are in the complex case). That is, given any open cover on a smooth manifold, one can find a partition of unity subordinate to that cover - this is a very topological property. Using partitions of unity you can paste together local functions as desired.
The existence of smooth partitions of unity comes down to the existence of a smooth (but not analytic!) bump function on $[-1,1]$. Edit: you can find details and formulas on Wikipedia.
A related fact is that on a paracompact smooth manifold, the sheaf of real-valued functions is fine (nLab,wikipedia).
A: This is not really an answer, but rather a puzzling    example about the notion of smoothness.
Consider  the unit disk 
$$ D=\lbrace z\in\mathbb{C};\;\;|z|\leq 1\rbrace.$$
The cyclic group $C_n:=\mathbb{Z}/n\mathbb{Z}$, $n\geq 3$,   acts on $D$ by rotations of angle  $2\pi/n$ about the origin.    Consider the quotient $D/C_n$.      The functions on this quotient  can be identified with the  $C_n$-invariant functions on $D$. There are   several  rings of functions on $D/C_n$
$$ C^0(D)^{C_n} :=\mbox{ continuous, complex valued  $C_n$-invariant functions on $D$}$$
$$ C^\infty(D)^{C_n} :=\mbox{smooth, complex valued  $C_n$-invariant functions on $D$}$$
$$ \mathcal{H}(D)^{C_n}:=\mbox{holomorphic valued  $C_n$-invariant functions on $D$}$$
Here is the surprise.  The ring $C^0(D)^{C_n}$ is isomorphic to the ring $C^0(D)$ of continuous complex valued functions on $D$. The reason  is that the  spaces $D$ and $D/C_n$ are homeomorphic compact  spaces and thus their rings of continuous complex valued functions  are isomorphic.    
Now observe  that the ring   $\mathcal{H}(D)^{C_n}$ is also  isomorphic to the ring  $\mathcal{H}(D)$. Indeed,  a holomorphic function
$$ f(z)=\sum_{k\geq 0} a_k z^k $$
is $C_n$ invariant iff   $a_k=0$,   $\forall  k\not\equiv 0  \bmod  n$. Thus
$$  f(z)\in \mathcal{H}(D)^{C_n} \Longleftrightarrow  f(z)=\sum_{k\geq 0} a_{kn} z^{kn} $$
The map  
$$ \mathcal{H}(D)\ni f(z) \mapsto  f(z^n)\in \mathcal{H}(D)^{C_n} $$
is the sought for isomorphism. These two examples   show that we cannot distinguish between $D$ and $D/C_n$   topologically or holomorphically.  Surprisingly 
$$ C^\infty(D)^{C_n} \not\cong C^\infty(D)$$
This is not obvious  but not terribly hard to prove. The upshot of this last fact is that  smoothly the disk $D$ and the cone $D/C_n$ are different.     
The point of this simple example is that smoothness is a rather subtle concept.   More    subtle examples can be found in   Kolmogorov's work on  Hilbert's 13th problem.  Using probabilistic ideas he  gives a precise quantitative   meaning to the  fact that smooth functions are fewer than  continuous functions. 

MR0112032 (22 #2890)  Kolmogorov, A.
  N.; Tihomirov, V. M. ε-entropy and
  ε-capacity of sets in function spaces.
  (Russian)  Uspehi Mat. Nauk 14 1959
  no. 2 (86), 3–86.
Vituškin, A. G.; Henkin, G. M.
  Linear superpositions of
  functions.  Uspehi Mat. Nauk 22
  1967 no. 1 (133), 77–124.

A: 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.
