function[F,Phi,Phi_Phi,Phi_F,Phi_R,BarRs,A_e, B_e,C_e]=mpcgain2(Ap,Bp,Cp,Nc,Np)
[m1,n1]=size(Cp);
[n1,n_in]=size(Bp);
%
A_e=eye(n1+m1,n1+m1);
%Forming the augmented model                          
A_e(1:n1,1:n1)=Ap;
A_e(n1+1:n1+m1,1:n1)=Cp*Ap;

B_e=zeros(n1+m1,n_in);
B_e(1:n1,:)=Bp;
B_e(n1+1:n1+m1,:)=Cp*Bp;

C_e=zeros(m1,n1+m1);
C_e(:,n1+1:n1+m1)=eye(m1,m1);

%Compute F matrix
h(1:m1,:)=C_e;
F(1:m1,:)=C_e*A_e;
for kk=2:Np
h((kk-1)*m1+1:kk*m1,:)=h((kk-2)*m1+1:(kk-1)*m1,:)*A_e;
F((kk-1)*m1+1:kk*m1,:)= F((kk-2)*m1+1:(kk-1)*m1,:)*A_e;
end

%Compute Phi matrix
v=h*B_e;
Phi=zeros(m1*Np,n_in*Nc); %declare the dimension of Phi
Phi(1:(m1*Np),1:n_in)=v; % first column of Phi

for i=2:Nc
Phi(:,((i-1)*n_in+1):(i*n_in))=[zeros((i-1)*m1,n_in);v(1:(Np-(i-1))*m1,:)];%Toeplitz matrix
end

BarRs=zeros(m1*Np,m1); 

if m1 == 1
    BarRs(((m1-1)*Np+1):m1*Np,m1)=1; % BarRs for SISO case [just vector of ones]
    else 
    for i= 2 : m1*Np/2+1 % BarRs for MIMO case [equivalent to Np blocks of eye(m1,m1) ]
    BarRs(2*i-3,1)=1;
    BarRs(2*i-2,2)=1;
    end 
end

Phi_Phi= Phi'*Phi;
Phi_F= Phi'*F;
Phi_R=Phi'*BarRs;