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.
It is a closed region, so max and min must occur. They can only occur on the boundary or at critical points of the function. So you can use the following steps:
Step 1: Find all the critical points of the function, and check whether they are in the constraint region.
Step 2: Use regular Lagrange multiplier method on the boundary of the disk.
Then combine the results from the two steps to find the max and min.
Best Answer
Because of the peculiar character of this function, we can make the substitution $ \ u \ = \ x + y \ $ to reduce the problem to finding the extrema of a function of a single variable, $ \ u^2 \ - \ u \ + \ 1 \ $. The domain to be studied becomes $ \ 0 \ \le \ u \ \le \ 2 \ $ . In the "interior" of the interval, we find, as godonichia remarks, that $ \ 2u \ - \ 1 \ = \ 0 \ \ \Rightarrow \ \ u \ = \frac{1}{2} \ $ , which corresponds to the extremal value $ \ \frac{3}{4} \ $ for the function. Upon investigating the endpoints of the interval, we find that the function takes on the value $ \ 3 \ $ at the endpoint $ \ u \ = \ 2 \ $ . This indicates that choice (c) is correct in the imaged question above.
The virtue of this approach is that constant values of $ \ u \ $ correspond to oblique lines of slope $ \ -1 \ $ passing through the first-quadrant unit square. This allows us to examine both the interior and the boundary of the square simultaneously, whereas a more usual investigation using the $ \ x- $ and $ \ y-$ coordinates requires us to check the interior, the four edges, and the four vertices of the square individually. (In this case, the edges can be groups into two pairs with the same behavior, and two of the four vertices has the same value of $ \ x + y \ $ . )