This might help to get beta as a function of x:
X = (0*2000/1100:0.1:2000/1100+10).';
n = 1e4;
trials = 1e3;
beta = zeros(numel(X),1);
Pr_failure = zeros(numel(X),1);
for j=1:numel(X)
pr_failure = zeros(trials,1);
x = X(j);
for i = 1:trials
RAND = rand(n,2);
% parameters
Mean_R = 2000 ; %mean
CoV_R = 0.1;%coefficient of variation
Std_R = Mean_R * CoV_R;
% Monte Carlo Simulation
Ri = Mean_R + Std_R .* norminv(RAND(:,1));
Mean_Q = x * 1100; % Determine which variable, x, will let the beta = target beta
CoV_Q = 0.18;
Std_Q = Mean_Q * CoV_Q;
Qi = Mean_Q + Std_Q .* norminv(RAND(:,2)); % simuate the Q term
g = Ri - Qi; % limit state function g = R - Q
m = find(g < 0); % count the number of failure cases
f = length(m); % find the failure cases number in the g = R-Q vector
pr_failure(i) = f / n; % probability of failure
end
pr_failure = mean(pr_failure);
pr_failure = max(eps,pr_failure);
beta(j) = norminv(1 - pr_failure); %Calculate the reliability index beta
Pr_failure(j) = pr_failure;
end
figure(1)
plot(X,beta)
figure(2)
plot(X,Pr_failure)
