Distributing n elements randomnly in r bins, and counting how many bins end up wih 1 element TL,DR
I will distribute n elements (indistinguishable) in r bins at random (equal chance of each bin) and try to calculate the expected value of X1="number of bins with exactly 1 element in it" as a function of n (r is known). Getting the probability distribution for X1 is interesting but not vital.    
Context
I'm trying to find the best strategy for a board game I play with my friends. In it, you can draw n cards (you decide how many beforehand) and pay an amount for each card you want to draw. There are r different types of card, all with equal chances of appearing. After this, getting only one card of a type is penalized, while getting a group of 2 or more of the same kind is rewarded. Obviously, drawing dozens of cards assures you will not get many "singletons", but the cost becomes prohibitive. Knowing how many singletons vs groups i can expect if i draw n cards can help me optimize my plays (yes, I'm THAT guy that can't stop optimizing irrelevant things).
The model
Each card type is a "bin", and each card drawn is an element. There are r bins. N elements are distributed in the bins. Each element has to go into one and only one bin (once you show the card, you know what type it is). Each bin has the same chances of being selected. There is no limit on how many elements a bin can hold. Elements are indistinguishable, because the only thing that matter is how many of each type you get, regardless of the order.
What I want to do
I'm trying to plot (expected number of bins with 1 element), (expected number of bins with 2 or more elements) vs (number of elements distributed), with a 10% and 90% percentile lines if possible. The graph tools i know all require a explicit formula (a function of n) for the plot lines and allow no code.
I know how to make a simulation of the card draw, but I prefer to use it as a validation for the analytical formula.
I'm very familiar with excel/google sheets and the rest of my related analisys are on that format. I can learn another tool if necesary, but the distribution has to be expressed as a function of n and r.
What I've tried
Reading the wikipedia article about probablility distributions, i find that this one https://en.wikipedia.org/wiki/Multinomial_distribution describes exactly the experiment I'm doing, but the information it provides is not the one I look for. This is already beyond the statistics I've learned in college and I'm stuck.
I also tried modelling the problem as a system of differential equations of n variables X0, X1...Xn where Xi="number of bins with exactly i elements in it". When distributing an element, the probability that it falls in a bin with i elements is proportional to Xi, and the effect is that Xi decreases in 1 and X(i+1) increases in 1. It looks similar to this:
dX0/dn=-k·X0
dX1/dn= + k·X0 - k · X1
dX2/dn= + k·X1 - k · X2
...
dXn/dn= + k · X(n-1) 
What failed in my approaches
About the multinomial, I have to provide the specific combinations of number of elements in each bin to get the chances of the combination, but I want to "group" all the combinations that have say, 3 bins with 1 element, and combine the chances of all those combinations. This becomes unwieldly when n increases.
About the differenctial equations,  I'm turning a discrete variable (number of elements) into continuous. It has undesirable effects, for example distributing only 1 element, and finding that there are 0,000001 bins with 20 elements already. Also, I only know the values at n=0 (all the bins are empty, so X0=r, and 0 for the rest) so I struggle giving a value to k.
My questions
• Is any probability distribution that models my problem and provides the information I look for better than the multinomial?
• If not, is it possible to combine explicitly all the combinations of / into ?
• What other data point (apart from n=0) can be used to estimate k in the case of the differential equations.
Observations
• n is at least 1 but has no upper bound 
• The cost of drawing a card is constant so it doesn't matter
• The amount of physical cards is the game is high enough that drawing a card without replacing it doesn't change the probabilities of getting each card type in any meaningful way, so chances of getting each type remain constant.
• First post here!  I really wanted to make a structured question, but maybe I was too verbose, sorry about that. Will try to correct it in later post once I get the hang of this.
 A: It may  interest the reader that  this problem can be  solved by total
enumeration and  inclusion-exclusion, which requires  the manipulation
of certain generating functions. We begin with the closed form result
$$\frac{1}{r^n}
\sum_{k=0}^r k {r\choose k} {n\choose k} k!
\sum_{q=0}^{r-k} {r-k\choose q} (-1)^q {n-k\choose q} q! 
(r-k-q)^{n-k-q}.$$
What this says is that we classify all configurations according to the
number $k$ of bins that contain  one ball, which may range from $0$ to
$r$ (first  sum). As we want  to compute the expectation  we include a
factor $k$  at the  start. Next we  choose the  $k$ bins from  the $r$
available ones (term  ${r\choose k}$). Then we choose  the trials from
among  the  $n$ trials  where  those  chosen  bins were  filled  (term
${n\choose  k} k!$).  We  then use  inclusion-exclusion to  obtain the
exact  count of  configurations corresponding  to the  choice  of bins
containing  one ball  by classifying  according to  the number  $q$ of
additional bins  containing one  ball.  We chose  these bins  from the
remaining  $r-k$  ones  (term   ${r-k\choose  q}$),  subtract  or  add
according to  inclusion-exclusion (factor $(-1)^q$)  and determine the
total count of configurations containing $k+q$ or more singleton bins,
where we choose the necessary $q$ trials from the $n-k$ remaining ones
(factor ${n-k\choose  q} q!$) and  distribute the rest  freely (factor
$(r-k-q)^{n-k-q}$).
Start the evaluation with the inner sum, getting
$$\sum_{q=0}^{r-k} {r-k\choose q} (-1)^q {n-k\choose q} q! 
(r-k-q)^{n-k-q}
\\ = (n-k)! \sum_{q=0}^{r-k} {r-k\choose q} (-1)^q 
\frac{(r-k-q)^{n-k-q}}{(n-k-q)!}.$$
Note however that
$$\frac{(r-k-q)^{n-k-q}}{(n-k-q)!}
= \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{\exp((r-k-q)z)}{z^{n-k-q+1}} \; dz$$
and we get for the sum
$$(n-k)! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{\exp((r-k)z) }{z^{n-k+1}} 
\sum_{q=0}^{r-k} {r-k\choose q} (-1)^q z^q \exp(-qz)
\; dz
\\ = (n-k)! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{\exp((r-k)z) }{z^{n-k+1}} 
(1-z\exp(-z))^{r-k} \; dz
\\ = (n-k)! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n-k+1}} 
(\exp(z)-z)^{r-k} \; dz.$$
On substituting this into the outer sum we get
$$\frac{n!}{r^n}
\sum_{k=1}^r k {r\choose k}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n-k+1}} 
(\exp(z)-z)^{r-k} \; dz
\\ = \frac{n!}{r^{n-1}}
\sum_{k=1}^r {r-1\choose k-1}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n-k+1}} 
(\exp(z)-z)^{r-k} \; dz
\\ = \frac{n!}{r^{n-1}}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}}
 (\exp(z)-z)^{r} 
\sum_{k=1}^r {r-1\choose k-1} \frac{z^k}{(\exp(z)-z)^k}
\; dz
\\ = \frac{n!}{r^{n-1}}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}}
 (\exp(z)-z)^{r-1} 
\sum_{k=0}^{r-1} {r-1\choose k} \frac{z^k}{(\exp(z)-z)^k}
\; dz
\\ = \frac{n!}{r^{n-1}}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}}
 (\exp(z)-z)^{r-1} 
\left(1+\frac{z}{\exp(z)-z}\right)^{r-1}
\; dz
\\ = \frac{n!}{r^{n-1}}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}}
\left(\exp(z)-z+z\right)^{r-1}
\; dz
\\ = \frac{n!}{r^{n-1}}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}} \exp(z(r-1))
\; dz.$$
Extracting the coefficient now yields
$$\frac{n!}{r^{n-1}} \frac{(r-1)^{n-1}}{(n-1)!}$$
which is
$$\bbox[5px,border:2px solid #00A000]{
n\left(1-\frac{1}{r}\right)^{n-1}}$$
