Problem in Monte Carlo simulation.

My code is not giving me the desired result. It always shows 1. For your ease let me explain the base of code again

The failure occur when "C >= Cr". so I intend to calculate the probabilty of that failure, lets say its P_fail. Here, C will be Calculated using

C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));

and Cr = 1.2

and ultimately final output is calculated when P_fail becomes greater than P_fail_Max (which equals to 0.10).

clear all
clc
n = 200;
%cover depth
mean_Cd = 16.1;
sigma_Cd = 0.23;
%diffusion coefficient at 28 days
mean_D = 3.87*10^(-12);
sigma_D = 0.52;
%surface chloride concn
mean_Cs = 13.1;
sigma_Cs = 0.103;
%Time exponent
mean_T = 0.2;
sigma_T = 0.2;
%critical chloride content
mean_Cr = 1.2;
sigma_Cr = 0.2;
%Now use normrnd function and lognrnd to generate values randomly
Cd = normrnd(mean_Cd, sigma_Cd, [n,1]);
D = lognrnd(mean_D, sigma_D, [n,1]);
Cs = normrnd(mean_Cs, sigma_Cs, [n,1]);
T = normrnd(mean_T, sigma_T, [n,1]);
Cr = normrnd(mean_Cr, sigma_Cr, [n,1]);
%Create for loop
for i=1:n
    rng(0,'twister'); %generating a vector of 200 random values for Cr values  
    a = 0.52;
    b = 3.87*10^(-12);
    D = a.*randn(200,1) + b;
    %display(D)
    rng(0,'twister'); %generating a vector of 200 random values for Cr values  
    a = 0.2;
    b = 0.2;
    T = a.*randn(200,1) + b;
    % display(T)
    if (D(i,1).*T(i,1) > 0 )  % Check to see if this number is real and positive before calculating sqrt
        con = D(i,1).*T(i,1);
    else
        con = 1;
        %break;
    end
    rng(0,'twister'); %generating a vector of 200 random values for Cr values  
    a = 0.2;
    b = 1.2;
    Cr = a.*randn(200,1) + b;
    %display(Cr)
    rng(0,'twister'); %generating a vector of 200 random values for Cr values  
    a = 0.23;
    b = 16.1;
    Cd = a.*randn(200,1) + b;
    % display(Cd)
    rng(0,'twister'); %generating a vector of 1000 random values for Cr values  
    a = 0.103;
    b = 13.1;
    Cs = a.*randn(200,1) + b;
    %display(Cs)
    % MAIN EQUATION TO CALCULATE C
    C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
    a = 1;
    b = i;
    C = a.*randn(20,1) + b;
    display(C)
    Pf = mean(C > Cr);
    display(Pf)
    
end
Pfmax = 0.10;
Tser = mean(Pf >= Pfmax);
display(Tser) % final output 

[abdellah Chougrani]

check the units

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/1707460-how-to-code-monte-carlo-simulation#comment_2181870