x1=-2:0.1:8;
x2=x1;
[X1,X2]=meshgrid(x1,x2);
S=4*X1-X2.^2-12;
H=25-X1.^2-X2.^2;
G1=X1.^2-10*X1+X2.^2-10*X2+34;
G2=-(X1-3).^2-(X2-1).^2;
mesh(X1,X2,S); hold;

for i=-20:20
 contour3(X1,X2,H+i,[i i]);
 contour3(X1,X2,G1+i,[i i]);
 contour3(X1,X2,G2+i,[i i]);
end

xlabel('x1'); ylabel('x2'); zlabel('S(x)');
axis([-2 8 -2 8 -20 20]);
view(20,50);

figure;

contour(X1,X2,S,20); hold
contour(X1,X2,H,[0 0]);
contour(X1,X2,G1,[0 0]);
contour(X1,X2,G2,[0 0]);
