% doa_analysis.m

% 2011/02/03   H. Ichihara

clear all;

% system
x = sdpvar(2,1);
f = [-x(2);x(1)+(x(1)^2-1)*x(2)];

% candidate Lyapunov function
v = 1.5*x(1)^2-x(1)*x(2)+x(2)^2;

% gamma
gam = 1.0; % 2.3044;

% epsilon
ep = 1.0e-6;

% SOS polynomial
z1 = [x(1);x(2)];
%s1 = qf_polynomial(z1);
s1 = polynomial(x,2,2);

% SOS condition
co = -jacobian(v,x)*f-s1*(gam-v)-ep*x'*x;

% SOS set
Fs = sos(s1)+sos(co);
params = recover(setdiff(depends(Fs),depends(x)));

% call solver
options = sdpsettings;
options.solver = 'sedumi'; % 'lmilab','sdpt3'
if strcmp(options.solver,'lmilab')
 % If you use LMILAB, you may need to change 'options.sos.model'.
 options.sos.model = 2;
end
[sol,monos,Q,resi] = solvesos(Fs,[],options,params)

% check set
for resii = resi'
 if resii > 1e-6
  disp('Infeasible!!')
 end
end

% EOF