th = sdpvar(2,1);
A = [-2.7-0.5*th(1) 1 -1; % system
     -12-2*th(2) 2 -1;
     3 1 -18.9-3.5*th(1)-2.7*th(2)-0.5*th(1)*th(2)];
n = length(A);
alpha = 1.0; ep = 1e-6;
g1 = alpha^2-th(1)^2; % constraints
g2 = alpha^2-th(2)^2;
zp = monolist(th,1);
P = sm_polynomial(zp,n); % Lyapunov variable
za = monolist(th,1);
zb = monolist(th,1);
Sa1 = qf_polynomial(za,n); % matrix SOS polynomials
Sa2 = qf_polynomial(za,n);
Sb1 = qf_polynomial(zb,n);
Sb2 = qf_polynomial(zb,n);
Ca = P-ep*eye(n); % SOS conditions
Ca = Ca-Sa1*g1-Sa2*g2;
Cb = -P*A-A'*P-ep*eye(n);
Cb = Cb-Sb1*g1-Sb2*g2;
Fs = sos(Sa1)+sos(Sa2)+sos(Sb1)+sos(Sb2)+sos(Ca)+sos(Cb);
params = recover(setdiff(depends(Fs),depends(th)));
[sol,monos,Q] = solvesos(Fs,[],[],params)