Solution (source code)

= Solution

Take $a,b,c>0$ as the semiaxis lengths. For a centered axis-aligned box with positive half-edge lengths $x,y,z$, the <volume> is $V=8xyz$. Maximize $xyz$ under $g=x^2/a^2+y^2/b^2+z^2/c^2=1$. The <Lagrange multiplier> equations are
$$
yz=2\lambda x/a^2,\qquad xz=2\lambda y/b^2,\qquad xy=2\lambda z/c^2.
$$
Multiplying by $x,y,z$, respectively, shows that $x^2/a^2=y^2/b^2=z^2/c^2$. The constraint makes each value $1/3$. All boundary boxes with a zero edge have zero volume, while the positive stationary box exists. <Compactness> therefore gives the <maximum box volume in an ellipsoid>
$$
\boxed{x=a/\sqrt3,\quad y=b/\sqrt3,\quad z=c/\sqrt3,\qquad V_{\max}=\frac{8abc}{3\sqrt3}.}
$$
Allowing a displaced or differently oriented rectangular box cannot improve this result. Let $A=\operatorname{diag}(a,b,c)$, let $c_0$ be its center, and let $u_1,u_2,u_3$ be its perpendicular half-edge vectors. Averaging $|A^{-1}(c_0+\sum\varepsilon_i u_i)|^2\le1$ over the eight signs gives $|A^{-1}c_0|^2+\sum|A^{-1}u_i|^2\le1$. The determinant bound $|\det(A^{-1}u_1,A^{-1}u_2,A^{-1}u_3)|\le\prod|A^{-1}u_i|$, followed by the <arithmetic-geometric mean inequality>, makes its volume at most $8abc/(3\sqrt3)$. The box found above attains the bound.