how do i properly code a formula with multiple summations in matlab?

hello this is my first ever question here,

somewhat new to using Matlab

i'm trying to recreate a figure of :The Average Uplink Spectral Efficiency as a function of the number of Bse Station antennas M for different channel models. i'm struggling particularly with the NLOS (no line of sight) scenario, it's formula contains multiple summations and when applied in matlab, i think my lack of understanding when it comes to arrays and nested loops made me unable to replicate the results desired, which are shown in the pics blow:

(the parameters needed are provided below the figure)

i feel such a formula should be rudimentry in application yet here i am struggling with it ...

FOR MORE CONTEXT, THE PAPER THIS IS FROM IS: (( Emil Björnson, Jakob Hoydis and Luca Sanguinetti (2017), “Massive

MIMO Networks: Spectral, Energy, and Hardware Efficiency”, Foundations and Trends R

in Signal Processing: Vol. 11, No. 3-4, pp 154–655. DOI: 10.1561/2000000093. ))

had no issues regarding the LOS and lower_bound plots so you can ignore those

below is my best attempt at ploting NLOS:

i apologize in advance for the bad grammer, and of course all help is appreciated greatly!

clear all
 
% parameters:
M = 1:100; % Number of BS antennas
 
SNR0_dB = 0; % SNR of the desired UE (in dB)
SNR0 = 10.^(SNR0_dB/10); % SNR of desired UE (in linear scale)
 
beta_bar_dB = -10; % Inter-cell interference strength (in dB)
beta_bar = 10.^(beta_bar_dB/10); % Inter-cell interference strength (in linear scale)

% % % Average SE in NLOS, formula (1.29) on page (184)

sum3=0;sum2=0;sum1=0;
% sum1= zeros(size(M));sum2= zeros(size(M));sum3= length(size(M));

    A = 1./( (1-1/beta_bar).^M );
    A = A-1;
      B = exp(1./(SNR0*beta_bar)).*expint(1./(SNR0*beta_bar));
        C = A.*B;
        C = C/log(2);
     
for m = 1:length(M)

     for L = 0:M-m
         
        for n = 1:L
            
            for j = 0:n-1
                sum3 = sum3 + 1/(factorial(j) .* (SNR0).^j) ;
                
            end
            sum2 = sum2 + (1/n)*sum3 ;
            sum3 = 0 ;
            
        end
        D = (exp(1./SNR0).*expint(1./SNR0));
        E = D + sum2;
        sum1 = sum1 +(( (-1).^(M-m-L+1) ).* E )./...
                  ( ((1 - 1./beta_bar).^m) .*...
                  ( factorial(M-m-L) .* SNR0.^(M-m-L) .* beta_bar.*log(2) ) );
                         
     end
     SE_NLOS =  C .* sum1;
    
end 

%  figure
plot(M, SE_NLOS,'green');

xlabel('Number of BS antennas (M)');
ylabel('Average SE (bits/s/Hz)');

[Torsten]

I don't understand how the graphics can be for SNR_0 = 0 dB while this variable appears in the denominator of various expressions.

[Fares Zaritt]

i remember i too was confused by this before. although (the signal to noise ratio) SNR0_dB = 0 dB, the SNR0 value used in calculation is in linear scale : SNR0 = 10^(SNR0_dB/10) = 1, we always convert from dB to linear scale ( when calculating using these formulas ), so the only way for SNR0 to equal 0 in calculation is for SNR0_dB to be a large negative value (e.g -∞)

in my code i did convert it at the start under parameters (also gave beta_bar the same treatment) :

SNR0_dB = 0; % SNR of the desired UE (in dB)

SNR0 = 10.^(SNR0_dB/10); % SNR of desired UE (in linear scale)

beta_bar_dB = -10; % Inter-cell interference strength (in dB)

beta_bar = 10.^(beta_bar_dB/10); % Inter-cell interference strength (in linear scale)

Answer 1

[Torsten]

% parameters:
M_array = 1:100; % Number of BS antennas
 
SNR0_dB = 0; % SNR of the desired UE (in dB)
SNR0 = 10.^(SNR0_dB/10); % SNR of desired UE (in linear scale)
 
beta_bar_dB = -10; % Inter-cell interference strength (in dB)
beta_bar = 10.^(beta_bar_dB/10); % Inter-cell interference strength (in linear scale)

% % % Average SE in NLOS, formula (1.29) on page (184)

sum3=0;sum2=0;sum1=0;
% sum1= zeros(size(M));sum2= zeros(size(M));sum3= length(size(M));
for kk = 1:numel(M_array)
    
M = M_array(kk);
    A = 1./( (1-1/beta_bar).^M );
    A = A-1;
      B = exp(1./(SNR0*beta_bar)).*expint(1./(SNR0*beta_bar));
        C = A.*B;
        C = C/log(2);
     
for m = 1:M

     for L = 0:M-m
         
        for n = 1:L
            
            for j = 0:n-1
                sum3 = sum3 + 1/(factorial(j) .* (SNR0).^j) ;
                
            end
            sum2 = sum2 + (1/n)*sum3 ;
            sum3 = 0 ;
            
        end
        D = (exp(1./SNR0).*expint(1./SNR0));
        E = D + sum2;
        sum2 = 0;
        sum1 = sum1 +(( (-1).^(M-m-L+1) ).* E )./...
                  ( ((1 - 1./beta_bar).^m) .*...
                  ( factorial(M-m-L) .* SNR0.^(M-m-L) .* beta_bar.*log(2) ) );
                         
     end

    
end 
SE_NLOS(kk) = C+sum1;
sum1 = 0;
end

%  figure
plot(M_array, SE_NLOS,'green');

xlabel('Number of BS antennas (M)');
ylabel('Average SE (bits/s/Hz)');

[Fares Zaritt]

thank you very much for your assistance with the code, i cant express my gratitute enough!