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 k^M words is equally probable.
Let x(k, M) 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:
<br />
x(k, M) = x(k - 1, M - 1) + x(k, M - 1)<br />
Obviously, x(k, M) = 0, \ k > M.
From here, we deduce that:
<br />
x(k, k) = x(k - 1, k - 1) = \mathrm{const}.<br />
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:
<br />
x(k, k + 1) = x(k - 1, k) + 1 \Rightarrow x(k, k + 1) = x(1, 2) + k - 1, \ k \ge 1<br />
Obviously, x(1, M) = 1, M \ge 1. This boundaty condition is consistent with the above recursion relation if and only if x(0, M) = 0, M \ge 0.
To reduce the recursion w.r.t. k, we perform a discrete Laplace transform w.r.t. k:
<br />
X_M(s) \equiv = \sum_{k = 0}^{\infty}{x(k, M) e^{-k \, s}}<br />
Multiply the recursion by e^{-(k - 1)s} and sum w.r.t. k from 1 to ∞, to get:
<br />
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<br />
<br />
X_{M + 1}(s) = (1 + e^{-s}) X_{M}(s), \ M \ge 1<br />
This recursion w.r.t. M is a geometric sequence with an easy solution:
<br />
X_M(s) = X_1(s) \, (1 + e^{-s})^{M - 1}, \ M \ge 1<br />
<br />
X_1(s) = \sum_{k = 0}^{\infty}{x(k, 1) \, e^{-k s}} = x(0, 1) + x(1, 1) e^{-s} = e^{-s}<br />
Then, take z \equiv e^{s}, and we have:
<br />
\tilde{X}_M(z) = X_M(s) = \frac{1}{z} \, \left(1 + \frac{1}{z} \right)^{M - 1}, \ M \ge 1<br />
Then, if you expand this binomial into a Laurent series near z = 0, you can identify the coefficients as x(k, M). We have:
<br />
\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<br />
So, I think the answer is:
<br />
P(k, M) = \left( \begin{array}{c} M - 1 \\ k - 1 \end{array}\right) \, \frac{1}{k^M}<br />