function Z = Kronecker(X,Y)

%Generate the symmetric Kronecker product of matrix X and Y.

[m,n]=size(X);
ntilde=n*(n+1)/2;

Z=zeros(ntilde);

for i=1:ntilde
    H=smat(unitvector(i,ntilde));
    Z(:,i)=0.5*svec(Y*H*X'+X*H*Y');
end