load build
 load tline

 A = full(A);
 n = length(A)

 for m=1:30
    B = randn(n,m);
    C = randn(m,n);

    [Us,A] = schur(A);
    B = Us'*B;
    C = C*Us; 

    U = lyapUR(A,B);

    id = [n:-1:1];
    L = lyapUR(A(id,id)',C(id)');
    L = L(id,id);

    P = U*U';
    Q = L*L';
   
    if m==1, figure(1), else, figure(2), end, clf
    svP = svd(P);
    svQ = svd(Q);
    Hsv = svd(P*Q);
    semilogy(svP,'b.','markersize',12), hold on
    semilogy(svQ,'r.','markersize',12)
    semilogy(Hsv,'k.','markersize',12)
    set(gca,'fontsize',16)
    legend('singular values of P','singular values of Q', 'Hankel singular values',1)
    xlabel('$k$','fontsize',16,'interpreter','latex')
    ylabel('$\sigma_k$','fontsize',16,'interpreter','latex') 
    title(sprintf('$m=p=%d$', m),'fontsize',16,'interpreter','latex')
    ylim([1e-10 5e5])
    pause
 end