Hello everyone i have this code, it generates 1 plot , now i want to have 10000 samples of it of different values by running a loop and then i need to find the pdf of it, how can i do that do that? please help me, thank you

clear; clc;

n_s = 1e4;

Y2S = @ (x) x * 365 * 24 * 60 * 60; % Transform year into seconds

T_year = 40; % in years

% Generate onset and formulation time

pd_t_form = makedist ('Normal', 'mu', 9.3, 'sigma', 3.162);

t_form_year = random (pd_t_form); % Generate crack formulation time: in year

t_form = Y2S (t_form_year); % in seconds

% Tube parameters

pd_t_p = makedist ('Normal', 'mu', 8.3, 'sigma', 1/3); % Pressure difference: Mpa

P_diff = random (pd_t_p);

pd_t_d = makedist ('Normal', 'mu', 22.23, 'sigma',. 5/3); % Outer diameter: mm

d = random (pd_t_d);

pd_t_t = makedist ('Normal', 'mu', 1.27, 'sigma',. 125 * 1.27 / 3); % inner diameter mm

thickness = random (pd_t_t);

% S_u = 713; % Ultimate tensile strength, MPa

% S_y = 362; % Yeild strength, MPa

T = Y2S (T_year);

% Scotte model parameters: the velocity unit is mm / s, for non-cold workded

% Alloy 600

alpha = 2.8e-9;

K_th = 9;

beta = 1.16;

F = .93;

a_0 = .1; % Initial crack length: mm

a_rup = 62.5; % threshold value mm

t_current = t_form;

sigma_stress = P_diff * d / 2 / thickness; % Stress factor

% Initial values ​​for recording the crack evolution process

a_k = a_0; % a_k: current crack length

t = t_current;

a = a_k;

time_step = 1;

index = 1; % Counter for saving the result

if t <T% If t_current does not exceed $ T $

while t

K = F * sigma_stress * sqrt (pi * a_k / 2); % Stress factor

a_kp = a_k + alpha * (K-K_th) ^ beta * time_step; % Calculate the crack length at the next time point

t_current = t_current + time_step; % Update the time point

a_k = a_kp; % Update current crack length

index = index + 1; % Update counter

time_step = time_step + 1;

if a_kp> a_rup% If the rupture happens

t (index) = t_current;

a (index) = a_kp;

break;

else if t_current> T

t (index) = T;

a (index) = a_kp;

break;

else

t (index) = t_current;

a (index) = a_kp;

end

end

end

else% If t_current exceeds $ T $

t = T; % Censored

end

end

plot (t, a)

xlabel ('Time'), ylabel ('Crack length')

Answer 1

[Ayush Gupta]

 The 1000 instances data of the variable can be stored in a global variable which can be later used for calculating a much better and approximate pdf. C is the global array and a is the variable for which we want to calculate the pdf of so after every iteration of the for loop we can add

C = [C a]; 

After getting C we can use fitdist function and fit a distribution accordingly. After this pdf function in MATLAB can used to find the pdf. Examples on how to use the pdf function can be accessed here.