Start with the equations that you have derived:
\begin{eqnarray*}
ye^{xy}&=&3\lambda x^2,\\
xe^{xy}&=&3\lambda y^2,\\
x^3+y^3&=&16.
\end{eqnarray*}
As a first step, show that none of $\lambda, x$ or $y$ can be zero. (If one of them is zero, then the first two equations show that all three must be zero, contradicting the third equation.) This is useful as we now know that we can divide by these terms at will.
Now comparing the first and second equations gives $x^3=y^3$, and since both are real, we get $x=y$. The third equation then gives $x=y=2$, and the first or second yields the value of $\lambda$. We find $f=e^4$ at this point, and it must be a maximum.
As regards the minimum, recall that the Lagrange multiplier method identifies the possible location of max/min points IF they exist. I don't think your example has a minimum: By taking $x^3 = N^3$ and $y^3=16-N^3$, where $N$ is a large positive number, the constraint $x^3+y^3=16$ is satisfied. But $xy\sim-N^2$ can be made arbitrarily large and negative, so that $e^{xy}$ can be made as close as we like to $0$, but of course would never equal zero. So we can say that the infimum
$$\inf\{e^{xy}:x^3+y^3=16\}=0,$$
but the minimum
$$\min\{e^{xy}:x^3+y^3=16\}$$
does not exist: for any candidate minimum at $(x_0,y_0)$, we can always find $(x_1,y_1)$ with $0<f(x_1,y_1)<f(x_0,y_0)$.
Another approach is to eliminate $y$ and to treat this as a single variable max/min problem for $h(x)=\exp[x(16-x^3)^{1/3}]$. This gives another way of understanding why there is no minimum point of the function.
![enter image description here](https://i.stack.imgur.com/T8Ucx.jpg)
There is already an accepted answer, but I thought I'd leave some remarks since this is sort of a curious constraint surface. The function $ \ f(x,y,z) \ = \ x^2 + y^2 + z^2 \ $ can of course be thought of as the squared-distance from the origin to a point on the surface $ \ x^3 + y^3 - z^3 \ = \ 3 \ $ . Since this is an "open" surface, the function $ \ f \ $ has no absolute maximum value. The discussion by user112018 and mathlove shows that the choice of coordinate value for the location of extrema then are either zero or given by
$$ \ \lambda \ = \ \frac{2}{3x} \ = \ \frac{2}{3y} \ = \ -\frac{2}{3z} \ \ \Rightarrow \ \ x \ = \ y \ = \ -z \ . $$
As mathlove has said, this offers eight possibilities, less one, since $ \ ( 0, 0, 0) \ $ is not on the surface. The three remaining categories are:
I -- two coordinates are zero
$$ -z^3 \ = \ 3 \ \ \Rightarrow \ \ (0, 0, -3^{1/3}) \ , \ y^3 \ = \ 3 \ \ \Rightarrow \ \ (0, 3^{1/3}, 0) \ , \ x^3 \ = \ 3 \ \ \Rightarrow \ \ (3^{1/3}, 0, 0) $$
$$ \Rightarrow \ \ f \ = \ 3^{2/3} $$
for all three points. These are found by following along each coordinate axis until it intersects the "dimple".
II -- one coordinate is zero
$$ y^3 \ - \ z^3 \ = \ 3 \ , \ y \ = \ -z \ \Rightarrow \ \ 2y^3 \ = \ 3 \ \ \Rightarrow \ \ (0, \left( \frac{3}{2} \right)^{1/3} , -\left( \frac{3}{2} \right)^{1/3} ) \ ,$$
$$ x^3 \ - \ z^3 \ = \ 3 \ , \ x \ = \ -z \ \ \Rightarrow \ \ ( \left( \frac{3}{2} \right)^{1/3} , 0 , -\left( \frac{3}{2} \right)^{1/3} ) \ ,$$
$$ x^3 \ + \ y^3 \ = \ 3 \ , \ x \ = \ y \ \Rightarrow \ \ 2x^3 \ = \ 3 \ \ \Rightarrow \ \ (\left( \frac{3}{2} \right)^{1/3} , \left( \frac{3}{2} \right)^{1/3} , 0 ) $$
$$ \Rightarrow \ \ f \ = \ 2 \cdot \left( \frac{3}{2} \right)^{2/3} \ = \ 2^{1/3} \cdot 3^{2/3} $$
for all three points. These are located along diagonals, in each of the three coordinate planes, meeting the "dimple" .
III -- no coordinate is zero
$$ x \ = \ y \ = \ -z \ \Rightarrow \ \ 3x^3 \ = \ 3 \ \ \Rightarrow \ \ (1, 1 , -1 ) \ \ \Rightarrow \ \ f \ = \ 3 \ . $$
This point lies on a "body diagonal" running to the "deepest corner of the dimple".
So, owing to the particular symmetries and anti-symmetries of this surface, the squared-distance function has three local and absolute minima [Category I] and what appear to be three additional (shallow) local minima [Category II] and a (very shallow) local maximum [Category III].
Best Answer
Have you tried the old-school approach? Every point $(x,y)$ on the unit circle can be written as $(\cos\theta,\sin\theta)$ for some $\theta\in[0,2\pi]$, and the stationary points of $\cos(\theta)\sin(\theta)^2$ lie at the roots of $\sin(\theta)\left(2\cos^2\theta-\sin^2\theta\right)$, hence the maximum is attained at $\theta=\arctan\sqrt{2}$ and its value is $\frac{2\sqrt{3}}{9}$. By symmetry ($(x,y)\mapsto(-x,y)$) the minimum is just the opposite.