MATLAB: How to simulate atom decay

MATLAB

I'm trying to use Matlab to simulate an atom decay process . I have only managed to simulate the single decay process.When it come to consecutive decay I got completely lost.
I attached my script for single decay below. Hope to get help here.
Thanks you.

Best Answer

N0 = 1e5;  % initial number of atoms
lambda1 = 0.0001; %decays/s
lambda2 = 0.00005;
dt = 100; % time step (in seconds)
M = 1500; % total number of time steps
t = 0:dt:(M-1)*dt; % time series
p1 = lambda1*dt;  % probability of having a decay in the time dt
p2 = lambda2*dt;
% Define the array that will be populated with the number of atoms at each
% time step
N1 = zeros(M,1); N2 = zeros(M,1); N3 = zeros(M, 1); 
N1(1) = N0;
tic
for j=2:M
    n_decaysN1toN2 = nnz(rand(N1(j-1),1) <= p1);
    n_decaysN2toN3 = nnz(rand(N2(j-1),1) <= p2);
      N1(j) = N1(j-1) - n_decaysN1toN2; % number of type 1 atoms that are left
      N2(j) = N2(j-1) + n_decaysN1toN2 - n_decaysN2toN3;
      N3(j) = N3(j-1) + n_decaysN2toN3;
end
toc
plot(t, [N1, N2, N3]), box off
legend({'# atoms type 1', '# atoms type 2', '#atoms type 3'})
legend boxoff
Related Question