Why does this code give error?

clear;clc;

fc = 3e8; 
Nb = 1000; %% Number of snapshots
c=3e8;
wavelength = 3e8/fc;     
d = 0.5*wavelength;  
N = 10; 
beta=2*pi/wavelength; % %Wavenumber             
A = 1;% %Signal amplitude
snr = 5;% %SNR(dB)
theta = [30 80];
sigma = sqrt((A^2)/(2*10^(snr/10)));% %Variance of noise
M = length(theta);% %Number of signals
a=(d)/(2*sin(pi/N));% %radius of the circular array

%source signal
for k=1:M 
    D(k,:) = randi(1,Nb);
    S(k,:) =A*(2*D(k,:) - 1);
end   

for i=1:M
    for k=1:N;
        phi(k)=(2*pi*(k-1))/N;
        SteeringVector(k,i)=exp(j*(beta*a*cos(theta(i)*pi/180)*cos(phi(k))-sin(theta(i)*pi/180)*cos(phi(k))));
    end
end
% White Gaussien noise
B = (sigma^2)*(randn(N,Nb)+j*randn(N,Nb))/sqrt(2);

%Array output:signal plus noise
X = SteeringVector*S+B;

% Estimation of the spatial correlation matrix of the observed signal
Rxx = X*X'/Nb;

theta1=[0:180];
for i=1:length(theta1)
 
    A1=zeros([N 1]);
    for k=1:N
       phiii(k)=(2*pi*(k-1))/N;
        A1(k,1)= exp(j*(beta*a*cos(theta1(i)*pi/180)*cos(phiii(k))-sin(theta1(i)*pi/180)*cos(phi(k))));
    end;
   PBeamforming (i)= real(diag(A1'*Rxx*A1))/(N^2);
end;


figure(1);
plot(theta1,10*log10(PBeamforming ));
title('Beamforming spectrum');
xlabel('Angle [degree]');
ylabel('PBeamforming [dB]');
grid on;

Answer 1

[the cyclist]

It's because

randi(1,Nb)

generates an Nb*Nb array. You need

randi(1,Nb,1)

as below

clear;clc;

fc = 3e8; 
Nb = 1000; %% Number of snapshots
c=3e8;
wavelength = 3e8/fc;     
d = 0.5*wavelength;  
N = 10; 
beta=2*pi/wavelength; % %Wavenumber             
A = 1;% %Signal amplitude
snr = 5;% %SNR(dB)
theta = [30 80];
sigma = sqrt((A^2)/(2*10^(snr/10)));% %Variance of noise
M = length(theta);% %Number of signals
a=(d)/(2*sin(pi/N));% %radius of the circular array

%source signal
for k=1:M 
    D(k,:) = randi(1,Nb,1);
    S(k,:) =A*(2*D(k,:) - 1);
end   

for i=1:M
    for k=1:N;
        phi(k)=(2*pi*(k-1))/N;
        SteeringVector(k,i)=exp(j*(beta*a*cos(theta(i)*pi/180)*cos(phi(k))-sin(theta(i)*pi/180)*cos(phi(k))));
    end
end
% White Gaussien noise
B = (sigma^2)*(randn(N,Nb)+j*randn(N,Nb))/sqrt(2);

%Array output:signal plus noise
X = SteeringVector*S+B;

% Estimation of the spatial correlation matrix of the observed signal
Rxx = X*X'/Nb;

theta1=[0:180];
for i=1:length(theta1)
 
    A1=zeros([N 1]);
    for k=1:N
       phiii(k)=(2*pi*(k-1))/N;
        A1(k,1)= exp(j*(beta*a*cos(theta1(i)*pi/180)*cos(phiii(k))-sin(theta1(i)*pi/180)*cos(phi(k))));
    end;
   PBeamforming (i)= real(diag(A1'*Rxx*A1))/(N^2);
end;


figure(1);
plot(theta1,10*log10(PBeamforming ));
title('Beamforming spectrum');
xlabel('Angle [degree]');
ylabel('PBeamforming [dB]');
grid on;

[Sadiq Akbar]

Thanks a lot dear @the cyclist.