Introducing the so called slack variables #s_1, s_2# the optimization problem is transformed into an equivalent one
Find local minima, maxima of
#f(x,y)=x^2y-xy#
subject to
#g_1(x,y,s_1) = x +y - s_1^2 - 3= 0#
#g_2(x.y,s_2)=x+ y + s_2^2-5=0#
The lagrangian is
#L(x,y,s_1,s_2,lambda_1,lambda_2) = f(x,y)+lambda_1 g_1(x,y,s_1)+lambda_2g_2(x,y,s_2)#
#L# is analytical so the stationary points include the relative maxima and minima.
The determination of stationary points is done solving for #x,y,s_1,lambda_1,s_2,lambda_2# the system of equations given by
#grad L(x,y,s_1,lambda_1,s_2,lambda_2) = vec 0#
or
#{(lambda_1 + lambda_2 - y + 2 x y=0),
(lambda_1 + lambda_2 - x + x^2=0),
( -3 - s_1^2 + x + y=0),
(-2 lambda_1 s_1=0),
(-5 + s_2^2 + x + y=0)
,(2 lambda_2 s_2=0)
:}
#
Solving we get
#((x = 0.451416, y = 2.54858, lambda_1 = 0.24764,
s_1 = 0., lambda_2 = 0., s_2 = 1.41421),
(x = 2.21525, y= 0.78475, lambda_1 = -2.69208, s_1= 0., lambda_2= 0., s_2= -1.41421),
(x = 0.472475, y= 4.52753, lambda_1= 0., s_1 = -1.41421, lambda_2 = 0.249242, s_2= 0.), (x = 3.52753, y = 1.47247, lambda_1= 0., s_1= -1.41421, lambda_2 = -8.91591, s_2 = 0.))#
The first and second points are at the boundaries of #g_1(x,y,0)=0# and the third and fourth at #g_2(x,y,0)=0# respectively.
Their qualification must be done with #f_(g_1)# and #f_{g_2}# respectively. So,
#f_{g_1}(x) = -(x-3) (x-1) x#
#f_{g_2}(x) =-(x-5) (x-1) x#
#d^2/(dx^2)f_{g_1}(0.451416)=5.2915# qualifying this point as a local minimum
#d^2/(dx^2)f_{g_1}(2.21525)=-5.2915# qualifying this point as a local maximum
#d^2/(dx^2)f_{g_2}(0.472475) =9.2915# qualifying this point as a local minimum
#d^2/(dx^2)f_{g_2}(3.52753) =-1.2915# qualifying this point as a local maximum
Attached a figure with a contour mapping with the points found.