How to explain the Maxwell Boltzmann distribution graph (physically)? While studying the kinetic theory of gases, I came to the section: "The Distribution of Molecular Speed" of a book. The book first (without any explanation) proposes the Maxwell's speed distribution law as
$$P(v) = 4\pi \left ( \frac {M}{2\pi RT} \right)^{3/2} v^2 e ^{-Mv^2/2RT} $$ and then shows the distributional graphs for various gasses at different temperatures as shown in the figure below.

One thing that bothered me was that it didn't provide any physical explanation for the shape of the curve. Therefore, I would like to know how and why, via molecular collisions, we achieve the same (exact) curve at a given temperature even if we start with a random distribution of molecular speed (but given that $ \bar K = \frac {3}{2}k_bT$)?
I would prefer an intuitive explanation rather than a mathematical one.

Note that you needn't to comment on this portion of the question if you wish. 
I searched for this over the internet and on Physics.SE but couldn't find any (do tell if you know of any). I also discussed this with my friend came with an ad-hoc explanation but that too raises many questions. If you are interested, the explanation is given below but as I mentioned, it is not a good one and has many loopholes.

We considered a distribution in which the temperature is $T$ and every molecule possesses the same speed (though in random direction) and hence the kinetic energy of each molecule is $\frac {3}{2}k_bT$. When we continue further we see that the molecules collide with one another and (most of the time) the collisions are between glancing and head-on collision. It turns out that some molecules collide with other molecules in such a way that in due time they gain kinetic energy, some other molecules collide in such a way that there is a net loss of their kinetic energy and hence this distribution is achieved.

But as I said earlier there are many questions to be asked and one them is given below.

  
*
  
*Why are some type of collisions favoured more over the others? 
  

 A: The phase space argument put forward by Shaswata is correct. From a thermodynamical perspective this is is merely the result of the of the three dimensional ideal gas relation
$$
P_{v} = \frac{n_{v}}{N_{total}} = \frac{1}{Z}
\exp\left({-\frac{E_{v}}{k_B T}}\right)
$$
where $E_{v}= \frac{m}{2}v^2$ is the energy,  $n_{v}$ is the number of particles with velocity $v=|\vec v|$, $N_{total} = \sum_v n_{v}$ is the total number of particles,  and 
$Z_ = \sum_v \exp\left({-\frac{E_{v}}{k_B T}}\right)$ 
is the partition function. You will find proper derivatives in many thermodynamics books. However, for me this is not an intuitive or conceptional argument. Thus, here is how I like to obtain the Maxwell-Boltzmann distribution. 


*

*Each component of the velocity vector, $\{v_x, v_y, v_z\}$, is normally distributed: This follows from the fact that we obtain a  normal distribution (= Gaussian), if we have "many" independent random processes which contribute "equally" to the observable variable.  For justification check out the Poisson theorem and how the Poissonian approached the normal distribution.  Or check out the so called Central limit theorem.

*The square of the velocity, $v^2 = \sum_{i=x,y,z} v_i^2$, is  distributed as $\chi^2_3$ (omitting the normalisation constant):
This simply is the definition of the $\chi^2$ distribution. Most often the  definition reads: If $Z_i \sim N(0, 1)$ then $Y=\sum_{i=1}^k Z_i^2 \sim \chi^2_{k}$. Note that the tilde "$\sim$" should be read as "is distributed as $\ldots$" and that the probability density function is given by $f_Y \propto \sqrt{Y} \exp\left(-Y/2\right)$, for $k=3$.

*However, we are not interested in distribution of 
$Y = \vec v^2/\sigma_v^2$, but we'd like to know the distribution of $v=|\vec v| = \sigma_v \sqrt{Y}$. Hence, we have to transform the $\chi_3^2$ distribution. However, this is straight forward. Using the transformation law (for continuous random variables)
\begin{align}
f_v \; dv &= f_y \cdot \left|\frac{dY}{dv}\right| \;dv
\\
&= f_y \cdot \frac{2 v}{\sigma_v^2} \; dv
\propto  v^2\exp\left( -\frac{v^2}{2\sigma_v^2}\right) \; dv
\end{align}
where $f_v$ is the so called probability density function for the random variable $v$,
we obtain the Maxwell Boltzmann distribution. 


In conclusion, the Maxwell Boltzmann distribution of $|\vec{v}|$ is a direct consequence of the fact that each velocity component, e.g. $v_x$, is normally distributed.
A: All the collisions are equally likely. The maxwell distribution can be derived from the Boltzmann probabilities- namely the probability of a molecule having energy $E$ is $e^{-\beta E}$. Once you can show this it is not so hard to derive the maxwell relation since for any velocity $v$ the number of states is proportional to $e^{-\beta \frac{1}{2}mv^2}$. It is also proportional to $4\pi v^2$ since in the 3D space $v_x,v_y,v_z$ can take any value given that $v_x^2+v_y^2+v_z^2=v^2$ which is basically proportional to the number of points on the surface of a sphere of radius $v$.
So the real question is how the particles achieve through constant molecular collisions the Boltzmann distribution. Why the probability of a molecule having energy $\frac{1}{2}mv^2$ is $e^{-\beta \frac{1}{2}mv^2}$.
Say in equilibrium the number of particles having velocity $v_i$ is $n_i$. We are restricted to a system where the $\sum n_i=N$ and the total energy $\sum n_i \frac{1}{2}mv_i^2=E$. The state having the maximum multiciplity is,
$$ \mathrm{arg\,max}_{n_1,n_2,\dots n_k}\frac{N!}{n_1!n_2!\dots n_k!} \qquad \sum n_i=N,\sum n_i \frac{1}{2}mv_i^2=E$$
We know that for the maximum multiplicity without the energy constraint is when all the $n_i$ are equal i.e. $n_i = N/k \quad\forall i$.
Hence it makes sense that when we have a very large energy ($v_j>v$) for one of the states the corresponding $n_j$ will be small. Otherwise, the other $n_i$ will be too small due to the energy constraint. And when we have one of the $n_i$ to be too large and the others are too small we have a small multiplicity.
Similarly, a low energy will correspond to a high number of molecules in that state.
Hence when you let all kinds of collisions to take place the state that is most likely ends up being the one in which most of the particles have low energy.
This can be seen from if you consider the collision of molecules as well. Say the two molecules have final velocities $v1,v_2$ and initial velocities as $u_1,u2$.
From the conservation of energy,
$$v_1^2+v_2^2=u_1^2+u_2^2=V^2$$
Basically, the solution space is the positive quarter of a circle of radius $V$. Now it is not hard to see that for $v_2>V/2$ the number of cases is proportional to $1/3$ (subtending an angle of $\pi/6$). On the other hand for $v_2<V/2$ we have 2 times more cases (subtending an angle of $\pi/3$). In other words $v_1^2$ and $v_2^2$ counterbalance each other. What happens as a result is that $v_1,v_2$ has a greater chance of reducing than of increasing.
A: Some of the answers provided above are pretty good. Let me try to provide an answer from a slightly different formulation for completeness. The approach I am going to take is to start from the general non-equilibrium description of the gas in question, and motivate what being in "equilibrium" means. The full non-equilibrium dynamics is governed by the Boltzmann equation. 
Physical Kinetics, the 10th book in the Landau-Lifshitz series is an excellent reference.
We start with a gas composed of particles obeying the dispersion relation $\epsilon = \frac{p^2}{2m}$ (the gas is isotropic) . In general, when the gas is not in equilibrium, it's physical state is described by the "Boltzmann distribution function" $f(x,p,t)$. Note that $\int{d^dp f(x,p,t)} = n(x,t)$, the number density at (x,t), d = number of spatial dimensions. We can see that this allows us to describe non-equilibrium situations through its dependence on $x$ and $t$. 
The distribution function described above is governed by the Boltzmann equation. Let me first write it down, and then I'll describe what it means:
$\scriptsize{ [\partial_{t}  + v.\nabla_{x}  + F.\nabla_{p}]f(x,p,t) = \\ \iiint d^dp_1 d^dp_2 d^dp_3 \delta^d(p+p_1 -p_2 - p_3) \delta(\epsilon(p) +\epsilon(p_1) - \epsilon(p_2) -\epsilon(p_3)) w(p,p_1 \to p_2,p_3) [f(p_2,x,t)f(p_3,x,t) - f(p,x,t)f(p_1,x,t)] }$.
First, the LHS. $v= \nabla_p (\epsilon)$ is the "velocity" of a particle with momentum $p$, whereas $F$ is the force experienced by it. The LHS is the net change in the number of particles in a given phase space volume due explicit variation on time (first term) , or, phase space flows out of the phase space volume (next two terms). The last two terms are usually called the "streaming terms".
Now, the "collison term": The RHS is the rate at which collisions contribute to changing $f$. There must be some collisions that result in a final state that lies in the given phase space volume. These contribute positively (+) to the collision term. In the reverse direction, there are collisions that involve a particle in the given phase space volume. These contribute negatively (-). We can motivate that the rate of any type of collision $q_A + q_B \to q_C + q_D$ must depend on the number of particles available to participate in this collision, ie, $f(q_A) f(q_B)$. These observations explain the term inside the square bracket.
A few more comments are in order. First, the delta functions. These come from the assumption that the dynamics of the gas is symmetric under space-time translations. This  implies that total energy and total momentum are conserved quantities, which can only be possible if we only allow, at the microscopic level, collisions that conserve energy and momentum. This explains the deltas.
Finally, the collision rate $w$. A natural question to ask is: why do we not have different collision rates for the forward and backward reactions? ie, why is $w(p,p_1 \to p_2,p_3) = w(p_2,p_3 \to p,p_1)$? The answer is that at the microscopic level, we assume that the dynamics is invariant under time reversal, hence, in particular, forward and backward reactions rates must be same. This is called the "principle of detailed balancing".
So, bringing the streaming terms on the RHS, we can see what the Boltzmann equation tells us; that the change in number of particles in a given phase space volume is totally accounted for by the rate at which particles are "streaming" in (or out), and the rate at which the particles are being "collided" into or out of, the given phase space volume.
Now, we are in a position to motivate what we mean by "equilibrium". Equilibrium is the situation when each of the above terms is explicitly zero. For the first term, this implies $\partial_t f=0$. For the streaming terms, this implies $\nabla_x f= 0$ and $F=0$. For the collision term, we see that this implies that $[f(p_2)f(p_3) - f(p)f(p_1)] = 0$. The $x,t$ dependence can be dropped, as now we are describing a gas in equilibrium.
We must remember the constraints imposed on the momenta due to the delta functions in the collision term. In general, if I want to impose the above condition, given the constraint $g(p)+g(p_1) = g(p_2)+g(p_3)$, for some function $g$, then one solution to the above equation (it can easily be checked) is $f(p) = e^{-\gamma g(p)}$, where $\gamma$ is some constant. Here, we have two things to conserve, energy and momentum, hence, the general equilibrium distribution is described by $f_{eq}(p) = e^{-\beta (\epsilon(p) - p.V)}$, where $\beta$ and $V$ are constants. Note that the case of non zero $V$ signifies a gas in uniform motion (a gas in uniform motion is also in "thermodynamic equilibrium") with velocity $V$. We will assume that we are in the rest frame of the gas at equilibrium, ie, $V=0$. Hence, we obtain $f_{eq} = e^{-\beta \epsilon(p)}$. We identify $T = \frac{1}{\beta}$ as the "Temperature" (Boltzmann constant $k_B$ has been set to 1)
To compute the equilibrium speed distribution, we must integrate over a shell of momentum magnitude $|p|$. The (d-1) dimensional area of said shell  ~ $ |p|^{d-1}$. Identifying $|p|$ as $m v$, where $m$ is the mass of a single particle, we finally obtain
$f_{Maxwell}(v)$ ~ $v^2 e^{-\beta \frac{p^2}{2m}}$ ~ $v^2 e^{- \frac{mv^2}{2T}}$, which is what we set out to show.
A: Semoi's answer is good. But since you say

I would prefer an intuitive explanation rather than a mathematical one

to express it more simply, and without formulae, the Maxwell distribution is the chi distribution with three degrees of freedom (the components of velocity in Euclidean space). The chi distribution is the distribution of the positive square root of the sum of squares of a set of independent random variables each following a standard normal distribution. The normal distribution is the result of the central limit theorem, which basically says that when you have loads of identical random variables added together, the result tends to a normal distribution.
So, Maxwell assumed that the underlying distribution of velocities should be random, and it should result from very large numbers of collisions, and consequently it should be normally distributed in each direction. The Maxwell distribution graph follows from that, using standard calculations.
Note on second part: I do not understand why you ask why some collisions are more favoured. All collisions just have the effect of randomising the velocity distribution both in magnitude and direction. They cannot change overall energy, which determines the mean of the normal distribution in each direction.
A: Spherical coordinates give R^2 sin
The difference between $$ dN/N \propto \exp(-v_x^2/s^2) dv_x $$ and $$dN/N \propto v^2 \exp(-v^2/c^2)dv$$ lies in the detail about correctly moving from $dv_x$ to $dv$. Consider the infinitesimal volume in cartesian and spherical coordinates $$dv_x dv_y dv_z \propto v^2 \sin\theta\, d\phi\, d\theta\, dv $$ This is how $v^2$ appears before exponent (plus to confirming $v_x^2+v_y^2+v_z^2=v^2$ inside it).
