% min S(x) = -6*x1 - 2*x2 + x1^2 + 2 * x2^2
%  s.a. g1(x) = 2*x1 + x2 - 8 <= 0
%       g2(x) = 4*x1 + x2 - 10 <= 0
%        x >= 0
%  x' = [1  6]

Q=[1 0;0 2]*2
c=[-6; -2]
A=[2 1;4 1]
b=[8; 10]
lb=[0; 0]
ub=[]

op=optimset('LargeScale','off');
[x,s,ex,out,lambda]=quadprog(Q,c,A,b,[],[],lb,ub,[],op)

[X1,X2]=meshgrid([0:0.1:4],[0:0.1:10]);
S=c(1)*X1+c(2)*X2 + (Q(1,1)*X1.^2 + Q(2,2)* X2.^2)/2;
l=contour(X1,X2,S,[-9.4 -9 -7 -5 -3 0 5 10 20 40 80 120]); hold on
clabel(l);
plot(X1(1,:),b(1)-A(1,1)*X1(1,:),'g',X1(1,:),b(2)-A(2,1)*X1(1,:),'r')
xlabel('x_1'); ylabel('x_2');
