Inner matrix dimensions must agree

I am writing a code in cvx using Matlab. But having a error :
Inner matrix dimensions must agree.

In expression

cvx_begin
variable W1(N+1,N+1) complex semidefinite;
expression R12;

R12=-rel_entr(t(2,1),t(2,1)+real (rho1 * trace
(si1* W1)))/log(2);


Where si1 is N+1*N+1 complex matrix.

And W1 as already defined variable W1(N+1,N+1) complex semidefinite, rho1 (constant)

But it is giving me a error :
Inner matrix dimensions must agree.

Do suggest me the modification.

function [rbar, t, tau, W1] =Optimize_R(T, eta, P1, N, rho1, rho2, rho3, si1, si2)
cvx_begin sdp
cvx_precision high
variable t(4) nonnegative;
variable tau(2) nonnegative;
variable W1(N+1,N+1) complex semidefinite; %W=t1*V1 >= 0
dual variables y1 y2 y3 y4;
expression R12;
expression R02;
expression R13;
expression R1;
expression R2;

%Achievable rate from WD1 to WD2 and The HAP in phase 2 are:
R12=-rel_entr(t(2,1),t(2,1)+real(rho1trace(si1W1)))/log(2); %equation 27
%xlog2(1+y/x)=-rel_entr(x,x+y)/log(2)
R02=-rel_entr(t(2,1),t(2,1)+real(rho2
trace(si1*W1)))/log(2); %equation 28

%Achievable rate of WD2 relaying info of WD1 to The HAP in phase 3:
R13=-rel_entr(t(3,1),t(3,1)+real(rho3tau(1,1)))/log(2); %equation 22
%Overall achievable Rates of WD1 and WD2 are:
R1=min(R12,R02+R13); %equation 12
R2=-rel_entr(t(4,1),t(4,1)+real(rho3
tau(2,1)))/log(2); %equation 23

rbar = min(R1, R2);                                                      %Auxilary Variable
maximize rbar;
subject to 
    t(1,1)+t(2,1)+t(3,1)+t(4,1)<=T;                                      %(1)
    t(1,1) >= 0;
    t(2,1)>= 0;
    t(3,1)>= 0;
    t(4,1)>= 0;
    tau(1,1)>= 0;
    tau(2,1) >= 0;
    diag(W1)==t(1,1);                                            %|[W]n,n|=t1 for n=1,2,...N+1  (26)
    W1 >= 0;  
    y1: -rbar<= R02+R13;
    y2: -rbar<=R12;
    y3: -rbar<=R2;  
    rbar>=0;
    y4: tau(1,1)+tau(2,1)<=eta*P1*real(trace(si2*W1));            %(29)
    W1==semidefinite(N+1);

cvx_end

Main Code:

%% Simulation parameters
        P1=1;           %HAP transmit power in t1 time, 30dBm=1W
        N0=10^-12 ;     %Noise power -90dBm
        eta=0.8;        %fixed energy harvesting effiency 0<eta<1
        T=1;            %Total time t1+t2+t31+t32
%Rayleigh fading Channel Coefficients
      beta=10^-2;                                    %path loss =[0 1] of order 10^-2 h~CN(0,Beta_K*I_N)
      a1=sqrt(beta/2)*(randn(1,1)+1i*randn(1,1));    %channel coefficient vector of HAP-WD1 1*1
      a2=sqrt(beta/2)*(randn(1,1)+1i*randn(1,1));    %channel coefficient vector of HAP-WD2 1*1 
      a12=sqrt(beta/2)*(randn(1,1)+1i*randn(1,1));   %channel coefficient vector of WD1-WD2 1*1  
      Nmax=50;
      pg=sqrt(beta/2)*(randn(3*Nmax,1)+1i*randn(3*Nmax,1));
      for N=5:5:50                                   %Number of Reflecting Surfaces
            g=pg(1:N).';                             %channel coefficient vector of HAP-IRS 1*N
            a1R=pg(Nmax+1:Nmax+N);                   %channel coefficient vector of IRS-WD1 N*1
            a2R=pg(2*Nmax+1:2*Nmax+N);               %channel coefficient vector of IRS-WD2 N*1  

%% Optimization for Problem P2 
% Optimal Solution for V2 using CVX
        gamma12=diag(a2R.')*a1R;
        [v2_bar]=Optimalv2_vector(N,gamma12,a12);
    % Optimal v2_bar = [v2,1 v2,2 ..... v2,N]
        Theta2_Optimal=diag(v2_bar);             %required phase shift matrix of IRS in phase-2

% Optimal Solution for V3 using CVX
        gamma2=diag(g)*a2R;
        [v3_bar]=Optimalv3_vector(N,gamma2,a2);
    % Optimal v3_bar = [v3,1 v3,2 ..... v3,N]
        Theta3_Optimal=diag(v3_bar);             %required phase shift matrix of IRS in phase-3

%% Optimization for Problem P4    
% Fixed parameters
    rho1=eta*P1*((abs(v2_bar*(diag(a2R.')*a1R)+a12))^2)/N0;
    si1=[diag(g)*a1R;a1]*[diag(g)*a1R;a1]';
    si2=[diag(g)*a2R;a2]*[diag(g)*a2R;a2]';
    rho2=eta*P1*(abs(a1))^2/N0;
    rho3=abs(v3_bar*(diag(g)*a2R)+a2)^2/N0;
    
% Optimizing problem P4 using CVX
    [R_bar,t_prim,tau_prim,W1_prim]=Optimize_R(N,T, eta, P1, rho1, rho2, rho3, si1, si2);
    V1_optimal=W1_prim/t_prim(1,1);

It will be less confusing it you stick with one user name.

You have a complicated mess, solving multiple problems, using for loop over N, etc. And as is now, it is not reproducible because it is missing data.

I suggest using a combination of whos for MATLAB variables and just typing CVX expressions at the command line (including MATLAB variable times a CVX variable) and look at the dimensions. Are they what they should be? Are the dimensions conformal? It appears something isn’t. Is rho1 actually a scalar?