Integer partition with fixed number of summands but without order For a fixed $n$ and $M$, I am interested in the number of unordered non-negative integer solutions to
$$\sum_{i = 1}^n a_i = M$$
Or, in other words, I am interested in the number of solutions with distinct numbers. For $n = 2$ and $M = 5$, I would consider solutions $(1,4)$ and $(4,1)$ equivalent, and choose the solution with $a_1 \ge a_2 \ge ... \ge a_n \ge 0$ as the representative of the class of equivalent solutions.
I know how to obtain the number of total, ordered, solutions with the "stars and bars" method. But unfortunately, I cannot just divide the result by $n!$ since that would only work if all the $a_i$ are distinct.

You want to know the number of partitions of $M$ into at most $n$ parts. A standard bijection (transposing the Young diagram) shows that this is equal to the number of partitions of $M$ into parts of size at most $n$. This number $p_n(M)$ has, for fixed $n$, generating function
$$\sum_{M \ge 0} p_n(M) t^M = \frac{1}{(1 - t)(1 - t^2)...(1 - t^n)}.$$
By computing the partial fraction decomposition of this rational function, you can write down a closed form for $p_n(M)$ (again, for fixed $n$). This is efficient in the regime where $M$ is large compared to $n$. I don't know what regime you care about. 
For what it's worth, the dominant term (for fixed $n$ as $M \to \infty$) is easy to extract: it's given by 
$$p_n(M) \approx \frac{1}{n!} {M+n-1 \choose n-1}$$
which follows from the fact that the dominant pole at $t = 1$ has multiplicity $n$ and from a computation of the coefficient of the corresponding term in the partial fraction decomposition. In other words, dividing the number you get from stars-and-bars by $n!$ is approximately correct (for fixed $n$ as $M \to \infty$) because in this regime the probability of any two of the numbers being the same becomes negligible. There is also a nice geometric way to see this, as $p_n(M)$ is just the number of non-negative integer solutions $x_2, ... x_n$ to
$$2x_2 + 3x_3 + ... + nx_n \le M$$
and this approximates the volume of the corresponding simplex in $\mathbb{R}^{n-1}$. 
See Wilf's generatingfunctionology for general background about generating functions and, for very powerful methods for extracting asymptotics, see Flajolet and Sedgewick's Analytic Combinatorics. 
