syms x y lambda real
f(x,y) = (x+y)^2;
g(x,y) = x^2 + y ^2 - 1;
L(x,y,lambda) = f - lambda*g
eqn1 = diff(L,x) == 0;
eqn2 = diff(L,y) == 0;
eqn3 = diff(L,lambda) == 0;
ss = solve([eqn1,eqn2,eqn3],[x,y,lambda]);
T = table(double(ss.x),double(ss.y),double(ss.lambda),double(f(ss.x,ss.y)));
T.Properties.VariableNames = {'x','y','lambda','f'}
fsurf(f,[-2 2 -2 2],'FaceAlpha',0.3);shg;view(0,90)
hold on
plot3(double(ss.x),double(ss.y),f(double(ss.x),double(ss.y)),'r*')
hold off
