
%Variable declaration
sdpvar x y;
Q=sdpvar(6);

%polynomial declaration
p=1+x^2-x*y^2+y^4;

%setting the vector of monomials
monomial_vector=[1 x y x^2 x*y y^2]';

%Setting the sos conditions as an sdp
Conditions=[Q>=0,coefficients(monomial_vector'*Q*monomial_vector-p,[x,y])==0];

%solve the sdp
solvesdp(Conditions,Q(3,3));

%display Q and its singular value decomposition
QQ=double(Q)
[U,S,V]=svd(QQ)

%ignore the last two rows of U since the corresponding eigenvalue is
%essentially zero
U=U*sqrt(S(:,1:4));

%reconstruct the factors
factors=monomial_vector'*U;
for i=1:4
    sdisplay(factors(i))
end

%check the sos
sdisplay(factors*factors');
