function ps22a_6242_2004(nv,nu)
% function ps22a_6242_2004(nv,nu)
%
% reserach function for problem 2.2a

m=max(nu,nv)+ceil(3*rand);  % dimension of the reduced system
n=2*m+30+ceil(30*rand);     % randomized number of states
A=randn(n);                 % generate A,B,C
B=randn(n,1);
C=randn(1,n);
s0=randn+j*randn;         % generate s0
Ai=inv(s0*eye(n)-A);      % (s0I-A)^{-1}

U0=zeros(m,n);            % to store U0=[C(s0I-A)^{-1};C(s0I-A)^{-2};...]
Ck=C;
for i=1:nu+1,                 
    Ck=Ck*Ai;
    U0(i,:)=Ck;
end
U0(nu+2:m,:)=randn(m-nu-1,n);

V0=zeros(n,m);       % to store V0=[(s0I-A)^{-1}B;(s0I-A)^{-2}B;...]
Bk=B;
for i=1:nv+1,                 
    Bk=Ai*Bk;
    V0(:,i)=Bk;
end
V0(:,nv+2:m)=randn(n,m-nv-1);
    
[U,S]=svd(U0',0);
[V,S]=svd(V0,0);
U=inv(U'*V)*U';
A1=U*A*V;
B1=U*B;
C1=C*V;
n1=size(A1,1);

N=2*(nu+nv+2);
e=zeros(1,N);
A1i=inv(s0*eye(n1)-A1);
Bk=B;
B1k=B1;
for i=1:N,
    Bk=Ai*Bk;
    B1k=A1i*B1k;
    y(i)=C*Bk;
    y1(i)=C1*B1k;
end
close(gcf)
bar(abs(y-y1)./(1+abs(y)));grid