function [X,Y,Z]=xyfun(N,A,B,Bd)

%% Model for example
alpha = 0;
% N = 8; 
% nfull = 16;
% nu = 2;
% [Gfull] = model_baltrunc(nfull);
% % make a discrete-time approximation
% Ts	= 1;
% Gfull = c2d(Gfull,Ts);
% [A B C Bd] = deal(Gfull.a,Gfull.b(:,1:nu),Gfull.c,Gfull.b(:,nu+1:end));
% %% Matrices dimensions (X,Y)
% % Matrices without 0 element
% X = cell((N),1);
% Y = cell((N),N);
% Z = cell((N),1);
%% Matrix X

for NN=1:N
    XN = 0;    
    for j=0:NN-1
        Aj = mpower(A,j);
        xn=(alpha^(NN-1-j))*Aj*Bd;
    % 	insert
        XN = XN +xn;
    end
    X{NN,1} = XN;
end
%% Matrix Y

% zeros in Y
for k=1:N
    for j=1:N
        Y{k,j}=zeros(size(B));
    end
end


for j=0:N-1
    Aj = mpower(A,j);
    yn=Aj*B;
% insert 1
    Y{j+1,1} = yn;
end

% insert 2
for i=1:N-1
    for j=1:N-1
        Y{i+1,j+1} = Y{i,j};
    end
end

%% Matrix Z

for j=1:N
    Ajj = mpower(A,j);
% 	insert
    Z{j,1} = Ajj;
end

% conversion
X = cell2mat(X);
Y = cell2mat(Y); 
Z = cell2mat(Z); 

