how to apply the monte carlo (MC) method on my BPSK_TX_RX code?

I made this code and I would like some help to apply the Monte Carlo (MC) method in this programming.

Help me!!!

the code follows below:

% BPSK
M = 2;                 

% número de bits ou símbolos
N = 1000000;   

% transmissor

% Gera bits aleatórios ( matriz de variaveis binarias)
Bits_ale = randi([0 M-1],1,N);

% Modulação BPSK 0 -> -1; 1 -> 1
Bits_bpsk = 2 * Bits_ale-1;        

% Gera bits de marca d'água aleatórios
Bit_wat = randi([0 M-1],1,N);

Theta = pi/8; 

b = 1/sqrt(2)*(randn(1,N) + 1i*randn(1,N)); % white gaussian noise with 0dB variance 

Eb_N0_dB = 0:20; % multiple Eb/N0 values

for k = 1:N 

    
    if Bit_wat(k)==1
        
        Bit_enviado(k) = Bits_bpsk(k) * exp(1i * Bits_bpsk(k) * Theta);
        
    else
        
        Bit_enviado(k) = Bits_bpsk(k) * exp(-1i * Bits_bpsk(k) * Theta);
    
    end
    
end


for iter = 1:length(Eb_N0_dB)
  
    % Noise addition
   
    d = Bit_enviado + 10^(-Eb_N0_dB(iter)/20)*b; % additive white guassian noise
  
    % receiver - hard decision decoding
   
    ipHat = real(d)>0;
   
    % counting the errors
   
    nErr(iter) = size(find(Bits_ale- ipHat),2);
    
end

simBer = nErr/N; % simulated ber

theoryBer = 0.5*erfc(sqrt(10.^(Eb_N0_dB/10))); % theoretical ber

% plot
close all
figure
semilogy(Eb_N0_dB,theoryBer,'b.-');
hold on
semilogy(Eb_N0_dB,simBer,'mx-');
axis([0 10 10^-5 1])
grid on
legend('theory', 'simulation');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for BPSK modulation');
hold on
grid on