Let $a_k(m,n)$ be the number of ways of placing $n$ distinct objects in $k$ distinct boxes if there must be at least $m$ objects in each box. Suppose first that $k=1$. Clearly $a_1(m,n)=0$ if $n<m$, and $a_1(m,n)=1$ if $m\ge n$, so the exponential generating function for the $a_1(m,n)$ is $$A_1(x)=\sum_{n\ge 0}a_1(m,n)\frac{x^n}{n!}=\sum_{n\ge m}\frac{x^n}{n!}=e^x-\sum_{n=0}^{m-1}\frac{x^n}{n!}.$$
Suppose now that you know $A_k(x)$ for some $k\ge 1$. To distribute $n$ objects amongst $k+1$ boxes you must first choose some of them to go into box $k+1$ and then distribute the remainder amongst the first $k$ boxes. Suppose that you put $r$ objects into box $k+1$. You can choose these $r$ objects in $\binom{n}r$ ways, and there are then $a_1(m,r)$ ways to ‘distribute’ them to box $k+1$ and $a_k(m,n-r)$ ways to distribute the remainder amongst the other $k$ boxes. Summing over the possible values of $r$, we find that $$a_{k+1}(m,n) = \sum_{r=0}^n\binom{n}r a_1(m,r)a_k(m,n-r)\;,$$ which is exactly what is needed to give us $A_{k+1}(x)=A_1(x)A_k(x)$. By induction, then, $$A_k(x)=A_1(x)^k=\left(e^x-\sum_{n=0}^{m-1}\frac{x^n}{n!}\right)^k\tag{1}.$$ If we set $$p_m(x)=\sum_{n=0}^{m-1}\frac{x^n}{n!},$$ the Taylor polynomial of order $m-1$ for the exponential function, $(1)$ can be written simply $$A_k(x) = \big(e^x-p_m(x)\big)^k.$$
Here is a very painstaking approach that may help you to see exactly what’s going on.
The possible values of $b_1$ are $0,1,2$, and $3$, so far starters we try
$$1+x+\frac{x^2}2+\frac{x^3}6$$
to account for $b_1$. Similarly, the possible values of $b_2$ are $1,2,3$, and $4$, so we try
$$y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}$$
to account for $b_2$. I’m using different indeterminates for now, because at this point I still need to keep the $b_1$ and $b_2$ contributions separate.
The product of these polynomials is
$$\begin{align*}
&y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}+\\
&xy+\frac{xy^2}2+\frac{xy^3}6+\frac{xy^4}{24}+\\
&\frac{x^2y}2+\frac{x^2y^2}4+\frac{x^2y^3}{12}+\frac{x^2y^4}{48}+\\
&\frac{x^3y}6+\frac{x^3y^2}{12}+\frac{x^3y^3}{36}+\frac{x^3y^4}{144}\;;
\end{align*}$$
however, we don’t want the terms in $x^ky\ell$ with $k\ge\ell$, since they correspond to having $b_1\ge b_2$. After we throw them away, we have
$$y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}+\frac{xy^2}2+\frac{xy^3}6+\frac{xy^4}{24}+\frac{x^2y^3}{12}+\frac{x^2y^4}{48}+\frac{x^3y^4}{144}\;.$$
Now replace $y$ by $x$, collect terms, and adjust the denominators to match the exponents to get
$$\frac{x}{1!}+\frac{x^2}{2!}+\frac{4x^3}{3!}+\frac{5x^4}{4!}+\frac{15x^5}{5!}+\frac{15x^6}{6!}+\frac{35x^7}{7!}\;,$$
which is the egf for boxes $1$ and $2$ combined. Multiply this by $e^{3x}$, and you’re done.
(And now that I’ve written this, I see that Markus has given you the abbreviated version of it.)
Best Answer
Consider an exponential generating function $F(x)$. The coefficient of $\frac{x^r}{r!}$ represents the number of ways to put some structure on $r$ labeled objects.
If $F$ and $G$ are two exponential generating functions, then the coefficient of $\frac{x^r}{r!}$ in $F(x) G(x)$ conveniently counts the number of ways to split $r$ labeled objects into two groups, then to put the structure counted by $F$ on the first group and structure counted by $G$ on the second group.
In this case, you need to split $r$ objects into a first group of things that are not put into any box, and a second group of things that are put into the $n$ boxes. You should find that $e^x$ counts the first structure and $\frac{1}{(1-x)^n}$ counts the second.
To see that the coefficient of $\frac{x^r}{r!}$ in $\frac{1}{(1-x)^n}$ counts the number of ways to put $r$ objects into $n$ boxes, with the order in each box mattering, you can apply the same trick. To put $r$ labeled objects in boxes in this way, you can first split the $r$ objects into $n$ groups, and then for each group count how many ways the objects can be put in the box. There are $k!$ ways to put $k$ objects into a single box, so the generating function counting this is $\sum_{k=0}^\infty k! \frac{x^k}{k!} = \frac{1}{1-x}$. Multiply this with itself $n$ times to get the function for $n$ boxes.