Label the urns from 1 to k. Then, each placement of the balls is equivalent to an M-letter word with the "letters" 1, 2, ..., k. Each of these [itex]k^M[/itex] words is equally probable.
Let [itex]x(k, M)[/itex] be the number of M letter words in which each of the k letters of the alphabet occurs at least once.
We can divide this set into two disjoint subsets: A set in which 1 enters exactly once and in which it occurs more than once. The first subset has x(k - 1, M - 1) elements. The second x(k, M - 1). Therefore, we have the recursion:
[tex]
x(k, M) = x(k - 1, M - 1) + x(k, M - 1)[/tex]
Obviously, [itex]x(k, M) = 0, \ k > M[/itex].
From here, we deduce that:
[tex]
x(k, k) = x(k - 1, k - 1) = \mathrm{const}.[/tex]
This constant is obviously equal to 1, since there is only one way of filling k urns with k balls (by placing exactly one ball in each urn).
Then, we have:
[tex]
x(k, k + 1) = x(k - 1, k) + 1 \Rightarrow x(k, k + 1) = x(1, 2) + k - 1, \ k \ge 1[/tex]
Obviously, [itex]x(1, M) = 1, M \ge 1[/itex]. This boundaty condition is consistent with the above recursion relation if and only if [itex]x(0, M) = 0, M \ge 0[/itex].
To reduce the recursion w.r.t. k, we perform a discrete Laplace transform w.r.t. k:
[tex]
X_M(s) \equiv = \sum_{k = 0}^{\infty}{x(k, M) e^{-k \, s}}[/tex]
Multiply the recursion by [itex]e^{-(k - 1)s}[/itex] and sum w.r.t. k from 1 to ∞, to get:
[tex]
e^{s} \, \left[ X_M(s) - x(0, M) \right] = X_{M - 1}(s) + e^{s} \, \left[ X_{M - 1}(s) - x(0, M - 1) \right], \ M > 1[/tex]
[tex]
X_{M + 1}(s) = (1 + e^{-s}) X_{M}(s), \ M \ge 1[/tex]
This recursion w.r.t. M is a geometric sequence with an easy solution:
[tex]
X_M(s) = X_1(s) \, (1 + e^{-s})^{M - 1}, \ M \ge 1[/tex]
[tex]
X_1(s) = \sum_{k = 0}^{\infty}{x(k, 1) \, e^{-k s}} = x(0, 1) + x(1, 1) e^{-s} = e^{-s}[/tex]
Then, take [itex]z \equiv e^{s}[/itex], and we have:
[tex]
\tilde{X}_M(z) = X_M(s) = \frac{1}{z} \, \left(1 + \frac{1}{z} \right)^{M - 1}, \ M \ge 1[/tex]
Then, if you expand this binomial into a Laurent series near [itex]z = 0[/itex], you can identify the coefficients as [itex]x(k, M)[/itex]. We have:
[tex]
\tilde{X}_M(z) = \frac{1}{z} \, \sum_{k = 0}^{M - 1}{\left( \begin{array}{c} M - 1 \\ k \end{array}\right) \, \frac{1}{z^k}} = \sum_{k = 1}^{M}{\left( \begin{array}{c} M - 1 \\ k - 1 \end{array}\right) \, \frac{1}{z^k}} \Rightarrow x(0, M) = 0, x(k, M) = \left( \begin{array}{c} M - 1 \\ k - 1 \end{array}\right), 1 \le k \le M, x(k, M) = 0, k > M[/tex]
So, I think the answer is:
[tex]
P(k, M) = \left( \begin{array}{c} M - 1 \\ k - 1 \end{array}\right) \, \frac{1}{k^M}[/tex]