You’re doing fine so far; you just haven’t gone far enough. First, while it’s true that you’re not interested in any of the terms $z^n$ with $n>r$, it’s easier to include them and look at $\left(\sum_{k\ge 0}z^{2k}\right)^8$, because then you can replace the infinite series by a simple generating function to get
$$\frac1{\left(1-z^2\right)^8}\;,$$
from which you want the coefficient of $z^r$. Note that this does no harm: the extra terms in the infinite series cannot contribute to the $z^r$ term anyway.
A useful generating function to know is
$$\frac1{(1-z)^{n+1}}=\sum_{k\ge 0}\binom{n+k}kz^k=\sum_{k\ge 0}\binom{n+k}nz^k\;.$$
Substitute $z^2$ for $z$ and take $n=7$, and you get
$$\frac1{\left(1-z^2\right)^8}=\sum_{k\ge 0}\binom{k+7}7z^{2k}\;.$$
In particular, the coefficient of $z^r$ is clearly $0$ if $r$ is odd, and if $r$ is even it’s $\dbinom{7+\frac{r}2}7$.
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
The generating function to ensure that the $i^\text{th}$ object appears at least $n+i$ times is as follows: \begin{equation} g(x) = \underset{1^{\text{st}} \text{ object}}{\underbrace{(x^{n+1}+x^{n+2}+\ldots+x^k)}}\underset{2^{\text{nd}} \text{ object}}{\underbrace{(x^{n+2}+x^{n+3}+\ldots+x^k)}}\ldots\underset{n^{\text{th}} \text{ object}}{\underbrace{(x^{2n}+x^{2n+1}+\ldots+x^k)}}. \end{equation}
Here, the power of $x$ in the first term of the product represents the number of times the first object is picked. Since the first object appears at least $n+1$ times, the smallest power of $x$ in the first term is $n+1$. The maximum number of objects to be chosen is $k$, and hence, the maximum power is $k$. Similarly, we get the later terms. Finally, the number of ways to choose $k$ objects is the coefficient of $x^k$ in the generator function $g(x)$.