// positive definite
L=diag([2 3]);
A=rand(2);A=(A+A')/2;
[V,d]=eig(A);
H=V*L*V';
x=-2:0.1:2;
[X1,X2]=meshgrid(x,x);
S=2+(H(1,1)*X1.^2+H(2,2)*X2.^2)/2+H(1,2)*X1.*X2;
figure; mesh(X1,X2,S); hold on;
xlabel('x1'); ylabel('x2'); zlabel('S(x)')
contour(X1(1,:),X2(:,1),S',2+[0.1 0.5 1 2 3 4 5 6 8 10 12]);

// negative definite
L=-diag([2 3]);
H=V*L*V';
S=14+(H(1,1)*X1.^2+H(2,2)*X2.^2)/2+H(1,2)*X1.*X2;
figure; mesh(X1,X2,S); hold on;
xlabel('x1'); ylabel('x2'); zlabel('S(x)')
contour(X1(1,:),X2(:,1),S',14-[0.1 0.5 1 2 3 4 5 6 8 10 12]);

// indefinite
L=diag([2 -3]);
H=V*L*V';
S=14+(H(1,1)*X1.^2+H(2,2)*X2.^2)/2+H(1,2)*X1.*X2;
figure; mesh(X1,X2,S); hold on;
xlabel('x1'); ylabel('x2'); zlabel('S(x)')
contour(X1(1,:),X2(:,1),S',14+[-10 -8 -6 -4 -2 -0.5 -0.1 0.1 0.5 2 4 6 8 10]);

// positive semi-definite
L=diag([2 0]);
H=V*L*V';
S=2+(H(1,1)*X1.^2+H(2,2)*X2.^2)/2+H(1,2)*X1.*X2;
figure; mesh(X1,X2,S); hold on;
xlabel('x1'); ylabel('x2'); zlabel('S(x)')
contour(X1(1,:),X2(:,1),S',2+[0.1 0.5 1 2 3 4 5 6 8 10 12]);
