A possible source of confusion here as to the use of the Lagrange multiplier $ \ \lambda \ $ appears to be due to the first Lagrange equation being factored incorrectly. The two equations are found properly, but we should have
$$ 3x^2 \ - \ 2 \lambda \ x \ = \ 0 \ \ \Rightarrow \ \ x \ ( \ 3x \ - \ 2 \lambda \ ) \ = \ 0 \ \ , $$
$$ 8y \ - \ 2 \lambda \ y \ = \ 0 \ \ \Rightarrow \ \ 2y \ ( \ 4 \ - \ \lambda \ ) \ = \ 0 \ \ . $$
So two of the possible solution sets are correctly given as $ \ ( 0 , \ \pm 1 ) \ $ and $ \ ( \pm 1, \ 0) \ $ by applying the constraint equation.
However, there is only the only value for the Lagrange multiplier found from the second equation, $ \ \lambda \ = \ 4 \ $ . This would be applied in the second factor of the first equation to obtain
$$ 3x \ - \ 2 \lambda \ = \ 0 \ \ \Rightarrow \ \ x \ = \ \frac{2 \cdot 4}{3} \ = \ \frac{8}{3} \ \ . $$
This raises an interesting issue in problems of this sort. In this case, the result produced by using the Lagrange multiplier value is inapplicable because there are no points with $ \ x \ = \ \frac{8}{3} \ $ on the constraint surface. [A similar circumstance arises in this problem.] Therefore, the extremal values of our function are
$$ f(0, \ \pm 1) \ = \ 4 \ \ [\text{absolute maximum}] \ \ ; $$
$$ f(1, \ 0) \ = \ 1 \ \ [\text{local maximum}] \ \ ; $$
$$ f(-1, \ 0) \ = \ -1 \ \ [\text{absolute minimum}] \ \ . $$
The "geometric" interpretation of this result is that the values for the function represent maximum and minimum values of the $ \ z-$ coordinate of the possible intersection points of the surface $ \ z \ = \ f(x,y) \ = \ x^3 \ + \ 4y^2 \ $ with the "vertical" cylinder $ \ x^2 \ + \ y^2 \ = \ 1 \ $ .
We show two views of the intersecting surfaces, with three of the four "critical points" visible; $ \ (0, \ -1, \ 4) \ $ is on the "far side" of the diagram at left.
$$ \ \ $$
The value found for the multiplier is the same regardless of the radius of that cylinder, which suggests that the coordinate $ \ x \ = \ \frac{8}{3} \ $ could become important if the constraint surface were "larger". If we were to use the cylinder $ \ x^2 \ + \ y^2 \ = \ 9 \ $ as the constraint "surface", for instance, our result from the Lagrange multiplier would be applicable, now giving us extremal values
$$ f(0, \ \pm 3) \ = \ 36 \ \ [\text{absolute maximum}] \ \ ; $$
$$ f(3, \ 0) \ = \ 27 \ \ [\text{local maximum}] \ \ ; $$
$$ f(-3, \ 0) \ = \ -27 \ \ [\text{absolute minimum}] \ \ ; $$
$$ f\left(\frac{8}{3}, \ \pm \frac{\sqrt{17}}{3}\right) \ = \ \frac{716}{27} \ \approx \ 26.52 \ \ [\text{local minima}] \ \ . $$
For these larger cylindrical radii, the intersection with the surface $ \ f(x,y) \ $ along the "positive $ -x \ $ " face of the cylinder develops a "ripple", leading to the appearance of further "critical points" with $ \ x \ = \ \frac{8}{3} \ $ .
The somewhat more complicated shape of the intersection of the surfaces is visible in the left-hand graph; the critical point $ \ (0, \ -3, \ 36) \ $ is not visible in these views.
The geometrical picture is the following: We are asked to find the local extrema of the distance from the point $(0,b)$ on the $y$-axis to points on the parabola $y=x^2$. From looking at a figure we can guess the following: If $b\gg1$ there are two local minima high up, and a local maximum at $(0,0)$. If $0<b\ll1$ there is just one local minimum at $(0,0)$, and the same holds when $b\leq0$.
The intended computation goes as follows: Set up the Lagrangian
$$\Phi:=x^2+(y-b)^2+\lambda(y-x^2)\ ,$$
and solve the system
$$\Phi_x=2x-2\lambda x=0,\quad \Phi_y=2(y-b)+\lambda=0,\quad y=x^2\ .$$
From $x(1-\lambda)=0$ we infer (i) $x=0$ or (ii) $\lambda=1$. In case (i) we then obtain $y=0$ and a certain value of $\lambda$, and in case (ii) we obtain $y=b-{1\over2}$. The condition $y=x^2$ then implies that case (ii) only leads to real solutions if $b\geq{1\over2}$, and in this case we have $x=\pm\sqrt{b-{1\over2}}$.
It follows that Lagrange's method has confirmed our geometric analysis of the problem. Note however that it is quite cumbersome to do a second derivative test in the framework of this method. Instead we can do the following: Consider the parametric representation $x\mapsto (x,x^2)$ of the parabola, and instead of $f$ plus constraint look at the pullback
$$\psi(x):=f(x,x^2)=x^2+(x^2-b)^2\qquad(-\infty< x<\infty)\ .$$
Now analyze this function $\psi$ as a function of one variable. You will get the same results (depending on $b$) as before, and in addition the second derivative test will confirm what you knew all along. The case $b={1\over2}$ is special: Here the first nonvanishing derivative is $\psi^{(4)}(0)=24$. Since $4$ is even and $24>0$ we have a local minimum there.
Best Answer
On two of the three boundaries, the function is just zero. On the last one, you know $x + y = 1$, so you can rewrite your function as a one-variable function
$$f(x, y) = f(x, 1 - x) = x^a e^{-x} (1 - x)^b e^{-(1 - x)} = \frac 1 e x^a (1 - x)^b$$
which isn't too hard to maximize.
Then for the interior, try to identify the critical points where $f_x = f_y = 0$.