function ps23cdef_6242_2004(n,m,M,R,K)
% function ps23cdef_6242_2004(n,m,M,R,K)
%
% Solves Problem 2.3, items (c)-(f)

if nargin<1, n=50; end
if nargin<2, m=10; end
if nargin<3, M=1; end
if nargin<4, R=0.2; end
if nargin<5, K=4; end

g=(2*(n^2)*K/M)*toeplitz([ 2 -1 zeros(1,2*n-1)]);
In=eye(2*n+1);
On=zeros(2*n+1);
b=[1;zeros(2*n,1)];
c=[zeros(1,n) 1 zeros(1,n)];
A=[On In; -g -(R/M)*In];
B=[zeros(2*n+1,1);b];
C=[c zeros(1,2*n+1)];
P=(M/2)*[g On;On In];

V=zeros(4*n+2,m);
Ai=inv(A);
Bk=Ai*B;
for i=1:m,
    Bk=Bk/norm(Bk);
    V(:,i)=Bk;
    Bk=Ai*Bk;
    Bk=Bk-V(:,1:i)*(V(:,1:i)'*Bk);
end
    
U=inv(V'*P*V)*V'*P;
A1=U*A*V;
B1=U*B;
C1=C*V;

A2=V'*A*V;
B2=V'*B;
C2=C*V;

w=linspace(0,100,1000);
g=squeeze(freqresp(ss(A,B,C,0),w));
g1=squeeze(freqresp(ss(A1,B1,C1,0),w));
g2=squeeze(freqresp(ss(A2,B2,C2,0),w));

close(gcf)
subplot(2,1,1);plot(w,real(g),w,real(g1)); grid
subplot(2,1,2);plot(w,real(g),w,real(g2)); grid