How find the maximum of the $x^3_{1}+x^3_{2}+x^3_{3}-x_{1}x_{2}x_{3}$ Let $$0\le x_{i}\le i,\, i=1,2,3$$ be real numbers. Find the maximum of the expression
$$x^3_{1}+x^3_{2}+x^3_{3}-x_{1}x_{2}x_{3}$$

My idea: I guess
$$x^3_{1}+x^3_{2}+x^3_{3}-x_{1}x_{2}x_{3}\le 0^3+2^3+3^3-0\cdot2\cdot 3=35$$ 
But I can't prove it. Can you help ?
This problem is a special case of:


Step 1:  rescale the variables (to see why this is valid, see the addition, below), replacing $y_1=x_1, y_2=x_2/2, y_3=x_3/3$ so that they are on $0\le y_i\le 1$.
on the interior we have
$$\nabla f= \langle 3y_1^2-y_2y_3, 3y_2^2-y_1y_3, 3y_3^2-y_1y_2\rangle=\langle 0,0,0\rangle$$
So

$$\begin{cases}3y_1^2=y_2y_3\quad (1) \\ 3y_2^2=y_1y_3\quad (2) \\ 3y_3^2=y_1y_2\quad (3)\end{cases}$$

Since the interior has $y_1y_2y_3\ne 0$ we can divide $(1)$ by $(2)$ to get ${y_1^2\over y_2^2}={y_2\over y_1}$ i.e. $y_1^3=y_2^3\iff y_1=y_2$ and by symmetry $y_1=y_2=y_3$, which means $2y_i^2=0$ impossible, hence no critical points on the interior.
The whole shape is a cube, so we have to check $6$ sides. On any of the planes where a $y_i=0$ we get a reduction to $y_j^3+y_k^3$ with $0\le y_i,y_j\le 1$ which you can solve by inspection to have the max value $2$.
The other possibilities are some $y_i=1$ on the boundary giving $y_j^3+y_k^3-y_jy_k+1$ with $0\le y_j, y_k\le 1$.
Optimizing this we need to solve
$$\langle 3y_j^2-y_k, 3y_k^2-y_j\rangle =\langle 0,0\rangle.$$
So

$$\begin{cases}3y_j^2=y_k \quad (1)'\\ 3y_k^2=y_j\quad (2)'\end{cases}$$

Again, dividing $(1)'$ by $(2)'$ we conclude $y_j=y_k=0$ which is not on the interior, so we discard it.
This just leaves the edge lines where $y_i=y_j=1$ and $0\le y_k\le 1$, this is a $1$-variable optimization problem:
We have $2+y_k^3-y_k$ i to be optimized, so that we check $3y_k^2-1=0\iff y_k=\sqrt{1\over 3}$, and the endpoint where $y_k=1$.
From here revert to $x$ coordinates and try the options out to get the max.

This method generalizes to any number of variables by induction, seeing the interior constraints will always reduce to all variables being $0$, so you can do the same sort of test after going back to normal coordinates. And if you do the computation even for just $2$ variables, you'll see how to compute without having to check them all.

Addition
For those wondering why it's sufficient to just consider the $y$ critical points instead of directly the $x$ ones, here is the justification:
Let
$$A=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3\end{pmatrix}$$
Then $f(y)=f(A^{-1}x)$ and $f(x)=f(Ay)$.
The critical points are where the derivative, $Df=0$ (in $x$ coordinates), that is to say:
$D_x(f(x)) = f'(Ay)$ by the chain rule since linear maps are their own derivatives. But since $A$ is invertible, $D_y(f(x))=f'(Ay)A=f'(x)A$. This is $0$ iff $f'(x)$ is, ergo we can test for critical points in either coordinate system.
It's the same idea why diffeomorphisms are Tangent bundle isomorphisms, you're doing the same process here with a simple, invertible linear map.
