function [Y,H2] = markovs(A,B,C,D,p)
% [Y,H2] = markovs(A,B,C,D,p)
% Form the matrix of p+1 Markov parameters , and compute it's H2 norm squared
% Y = [ D CB CAB CA^2B CA^3B ...  CA^(p-1)B ]
% for the discrete time system , 
%  x(k+1)  =  A x(k)   + B u(k) 
%  y(k+1)  =  C x(k+1) + D u(k+1)

 [m,r] = size(D);
 if nargin < 5, help(markovs); return; end
 Y = zeros(m, (p+1)*r);                       % dimension of Y
 Y(:,1:r) = D;                                % first block  is D
 AkB = B;
 Y(:,r+1:2*r) = C * AkB;                      % second block is C*B
 H2 = norm(D,'fro')^2 + norm(C*B,'fro')^2;    % H2^2 norm of [D] and [CB]
 for k=2:p-1
   AkB = A * AkB;
   Y(:,k*r+1:k*r+r) =  C * AkB;
   H2 = H2 + norm(C*AkB,'fro')^2;             % add in H2^2 norm of C A^{k-1} B
 end

% --------------------------------------------------------------------- MARKOVS
