Найти экстремум функции при наложеных ограничениях

1 способ (символьно)

z= x y;
Maximize[{z,x+y==1},{x,y}]

2 способ (численно)

z= x y;
NMaximize[{z,x+y==1},{x,y}]

Визуализация решения

Clear["Global`*"]
Quiet[Solve[{x+y==1,z==x y},{x,y,z}]]
Show[
Plot3D[x y,{x,-1,1},{y,-1,1},AspectRatio->1,MeshFunctions->{#1+#2&},MeshStyle->Gray],
ParametricPlot3D[{x,1-x,x-x^2},{x,-1.1,1.1},PlotStyle->Directive[Thickness[0.01],Black]],
Graphics3D[{Black,Sphere[{1/2,1/2,1/4},0.04]}]
]

Clear["Global`*"]
u=(x1 x2^2)^(1/3);
M=p1 x1+p2 x2;
Minimize[{u,M>0,p1>0,p2>0,x1>0,x2>0},{x1,x2}]

p1=1;
p2=1;
u=(x1 x2^2)^(1/3);
M=p1 x1+p2 x2;
Plot3D[{u,M},{x1,-1,10},{x2,-10,10}]