By way of enrichment here is an approach using generating
functions. The generating function for this problem is
$$f(u, x, y) = \left(u + x+y+xy\right)^{10}.$$
The variable $u$ marks empty boxes. We get for the total number of
outcomes
$$[x^5] [y^5] f(1, x, y) =
[x^5] [y^5] (1+x)^{10} (1+y)^{10}
\\ = {10\choose 5}^2 = 63504.$$
Now we are interested in the coefficients on $[u^0]$ (no empty box),
$[u^1]$ (one empty box) and $[u^2]$ (two empty boxes). We extract
these coefficients by writing
$$f(u, x, y) = \sum_{q=0}^{10} {10\choose q}
u^q \left(-1 + (1+x)(1+y)\right)^{10-q}.$$
This yields for $[u^0]$
$$[x^5][y^5] \sum_{p=0}^{10} {10\choose p} (-1)^{10-p}
(1+x)^p (1+y)^p
\\ = \sum_{p=0}^{10} {10\choose p} (-1)^{10-p} {p\choose 5}^2
\\ = 252.$$
This is of course ${10\choose 5}$ as with no empty box every box
contains one ball and hence it remains to choose the location of the
five balls from the first round. Continuing we get for $[u^1]$
$$10 [x^5][y^5] \sum_{p=0}^{9} {9\choose p} (-1)^{9-p}
(1+x)^p (1+y)^p
\\ = 10 \sum_{p=0}^{9} {9\choose p} (-1)^{9-p} {p\choose 5}^2
\\ = 6300.$$
Finally we have for $[u^2]$
$$45 [x^5][y^5] \sum_{p=0}^{8} {8\choose p} (-1)^{8-p}
(1+x)^p (1+y)^p
\\ = 45 \sum_{p=0}^{8} {8\choose p} (-1)^{8-p} {p\choose 5}^2
\\ = 25200.$$
We thus have for the result
$$\frac{252+6300+25200}{63504} = \frac{1}{2}.$$
Addendum. We can ask about the number of configurations with
$m$ empty boxes when we have two rounds of $n$ balls being distributed
into $2n$ boxes. We obtain for the generating function
$$f(u, x, y) = \left(u + x+y+xy\right)^{2n}$$
which we write as
$$f(u, x, y) = \sum_{q=0}^{2n} {2n\choose q}
u^q \left(-1 + (1+x)(1+y)\right)^{2n-q}.$$
We thus obtain for the coefficient on $[u^m]$
$${2n\choose m} [x^n][y^n]
\sum_{p=0}^{2n-m} {2n-m\choose p}
(-1)^{2n-m-p} (1+x)^p (1+y)^p
\\ = {2n\choose m}
\sum_{p=n}^{2n-m} {2n-m\choose p}
(-1)^{2n-m-p} {p\choose n}^2.$$
Now to simplify this note that
$${2n-m\choose p} {p\choose n}
= \frac{(2n-m)!}{(2n-m-p)! n! (p-n)!}
= {2n-m\choose n} {n-m\choose p-n}$$
and we thus obtain
$${2n\choose m} {2n-m\choose n}
\sum_{p=n}^{2n-m} {n-m\choose p-n}
(-1)^{2n-m-p} {p\choose n}
\\ = {2n\choose m} {2n-m\choose n}
\sum_{p=0}^{n-m} {n-m\choose p}
(-1)^{n-m-p} {p+n\choose n}.$$
Using the Egorychev method we have
$${p+n\choose n} =
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}} (1+z)^{p+n} \; dz$$
and we get for the sum
$${2n\choose m} {2n-m\choose n} (-1)^{n-m}
\\ \times \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}} (1+z)^{n}
\sum_{p=0}^{n-m} {n-m\choose p} (-1)^p (1+z)^p
\; dz
\\ = {2n\choose m} {2n-m\choose n} (-1)^{n-m}
\\ \times \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}} (1+z)^{n}
(1-(1+z))^{n-m}
\; dz
\\ = {2n\choose m} {2n-m\choose n}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}} (1+z)^{n}
z^{n-m}
\; dz
\\ = {2n\choose m} {2n-m\choose n}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} (1+z)^{n}
\; dz
\\ = {2n\choose m} {2n-m\choose n} {n\choose m}.$$
We can confirm these with a simple Maple script up to about $n=7$
which goes as follows:
with(combinat);
GFENUM :=
proc(n)
option remember;
local res, c1, c2, empty, filled;
res := 0;
for c1 in choose(2*n, n) do
for c2 in choose(2*n, n) do
filled := {op(c1), op(c2)};
empty := 2*n-nops(filled);
res := res + u^empty;
od;
od;
res;
end;
GFX :=
proc(n)
add(binomial(2*n,m)*binomial(2*n-m,n)*binomial(n,m)*u^m,
m=0..n);
end;
E.g. we get for $n=6$ the OGF
$$924\,{u}^{6}+33264\,{u}^{5}+207900\,{u}^{4}+369600\,{u}^{3}
\\+207900\,{u}^{2}+33264\,u+924.$$
