Optimize function on $x^2 + y^2 + z^2 \leq 1$ Optimize $f(x,y,z) = xyz + xy$ on $\mathbb{D} = \{ (x,y,z) \in \mathbb{R^3} : x,y,z \geq 0 \wedge x^2 + y^2 + z^2 \leq 1  \}$. The equation $\nabla f(x,y,z) = (0,0,0)$ yields $x = 0, y = 0, z \geq 0 $ and we can evaluate $f(0,0,z) = 0$.
Now studying the function on the boundary $x^2 + y^2 + z^2 = 1$ gets really hairy. I tried replacing $x$ with $\sqrt{1 - y^2 - z^2}$ in order to transform $f(x,y,z)$ into a two-variable function $g(y,z)$ and optimize it on $y^2 + z^2 \leq 1$ but $g(y,z)$ is a pain to differentiate. I then tried spherical coordinates which really did not make it any much easier.
Got any suggestions on how to tackle it?

Another way:
$f$ is convex in each variable $x, y, z$ so the maximum can be obtained only on the boundary.  This means $x^2+y^2+z^2=1$.  Further,
$$4 \,f = 4xy(1+z) \leqslant 2(x^2+y^2)(1+z) = 2(1-z^2)(1+z) = (1+z)(2-2z)(1+z)$$
The RHS now is a product of three positive terms with a constant sum, so it gets maximised when all terms are the same, viz. when $1+z = 2- 2z \implies z=\frac13$.  As the first inequality becomes an equality iff $x=y$, the maximum is indeed achieved when $z = \frac13, \; x=y=\frac23$ and $\displaystyle f_{max} = \frac{16}{27}$.
