Simulation for buffer energy model

I am trying to compare my montecarlo simulation result with the baseline model. I have the result and figure which I have got from the monte carlo simulation. And now I am trying to create a realistic energy model, calculate energy outage probability and compare the same.

I am not sure I am creating the energy model and then calculating enery outage correctly (last two line of the code). Can anyone please help me with this. Thank You!

clc
clear all;

%Parameters defined

NoB=1;                  % Noise Density
%rt=2*1.4426;           % Rate converted into bits/sec/Hz---(1 nat = 1/ln
                          %2 bits) and (1 bits = 1/log e nats)
rt=3.0;                   % Rate in nats/sec/Hz
EH=0;                   % Initialising Energy Harvesting
Bmax=50;                % Threshold level of battery in Joules

PmaxdB=[15];
Pmax=10.^(PmaxdB/10);

Qu=[-3:0.1:0];
QoS_Component_u=10.^(Qu);



for i=1:length(QoS_Component_u)
    u=QoS_Component_u(i);
        
  EHH(i)=0;                       %Initialising EHH equal to zero
    muo_sum(i)=0;                   %Initialising mu_o equal to zero
    
    for j=1:100000                  %Monte-Carlo Simulations
        

        hx=(1/sqrt(2))*randn(1,1);
        hy=(1/sqrt(2))*randn(1,1);
        ht_rayleigh=hx.^2+hy.^2;  
          
        muo(i)=min((((exp(rt)-1)*NoB)./ht_rayleigh),Pmax); 
        EHH(i)=EHH(i)+exp(u*muo(i));
       
    end

    EHHH(i)=(EHH(i))/100000;            %Average for Expectation
    mufix(i)=(1/u)*log(EHHH(i));  
   
    mu_o_avg(i)=(muo_sum(i))/100000;    %Average of mu_o(i)
    
    Epsilon(i)=(mu_o_avg(i))/(mufix(i));
   
    Pout(i)= (Epsilon(i)*(exp(-u*(Bmax-Pmax))));    % Calculating Pout
    
    energy_left(i) = (Bmax + mufix(i)) - muo(i);   % creating energy model
    el(i)= energy_left(i)/100000;                  % calculate energy outage probability
end

[Dhawal Beohar]

@Torsten I really appreciate if you can have a look.

Answer 1

[Yash]

There's missing logic in the Monte Carlo simulation. muo_sum is initialized but never updated.

muo_sum(i)=0; %Initialising mu_o equal to zero
...
mu_o_avg(i)=(muo_sum(i))/100000; %Average of mu_o(i)
...

mu_o_avg(i) will always be evaluated as 0 and hence subsequent calculations of Epsilon, Pout and other parameters will also be incorrect.