Solution (source code)

= Solution

Expand both sides through second order as $\epsilon\to0$, holding $X,Y$ fixed. The product is
$$
e^{\epsilon X}e^{\epsilon Y}=I+\epsilon(X+Y)+\epsilon^2\left(\frac{X^2}{2}+XY+\frac{Y^2}{2}\right)+O(\epsilon^3).
$$
Put $Z=\epsilon(X+Y)+\epsilon^2[X,Y]/2$. Its exponential is
$$
\begin{aligned}
e^Z&=I+\epsilon(X+Y)+\frac{\epsilon^2}{2}\{[X,Y]+(X+Y)^2\}+O(\epsilon^3)\\
&=I+\epsilon(X+Y)+\epsilon^2\left(\frac{X^2}{2}+XY+\frac{Y^2}{2}\right)+O(\epsilon^3).
\end{aligned}
$$
Thus
$$
\boxed{e^{\epsilon X}e^{\epsilon Y}=e^{\epsilon X+\epsilon Y+\epsilon^2[X,Y]/2}+O(\epsilon^3).}
$$
This is the quadratic part of the <Baker--Campbell--Hausdorff formula>. The remainder estimate is in any fixed finite-dimensional <matrix> <norm>, and is an asymptotic statement rather than an exact formula for arbitrary finite $\epsilon$.