# Use Monte Carlo method to simulate consecutive decay in MATLAB

I'm trying to use MATLAB to simulate an atom decay process by using Monte Carlo approach. The process is as follow:

Suppose that atom 1 decay to atom 2, which, in turn, decay to stable atoms of type 3. The decay constants of 1 and 2 are $\lambda_1$ and $\lambda_2$, respectively.

Assume that at $t = 0$, $N_1 = N_0$ and $N_2 = N_3 = 0$.

$\lambda_1=0.0001$, $\lambda_2=0.00005$, $t_{\rm MAX}$ for the process is $150000\rm s$.

I'm not really familiar with Monte Carlo method, so 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.

N0=input('Enter total number of nuclides: '); %Determine N0
lambda1 = 0.0001; %decays/s
dt = 100; %time step
M = 1500; % total number of time steps
t = 0:dt:(M-1)*dt; %time series
p = lambda1*dt; %probability of having a decay in the time dt

% Define the array that will be populated with the number of atoms at each time step N = zeros(M,1);

N(1) = N0;

for j=2:M % Loop over time
nd = 0; % number of decays
for i=1:N(j-1)
%generate a single random number
r = rand();
%binary value =1 if the condition is True -> Atom has decayed
rb = r < dt*lambda1;
if rb == 1
nd = nd+1;
end
end
N(j) = N(j-1)-nd; % number of atoms left at the step 2
end

plot(t, N)


The decay of radioisotopes can be described by a set of coupled first order differential equations $$\frac{dN_{i}(t)}{dt} = -\lambda_{i}N_{i}(t) + \sum_{j \neq i}\lambda_{j\to i} N_{j}(t)$$ where $\lambda_{j\to i}$ represents the decay constant associated with the decay of isotope $j$ to isotope $i$, $\lambda_{i}$ is the total decay constant of isotope $i$ and $N_{i}(t)$ is the time-dependent concentration of isotope $i$.
This can be easily written in matrix form $$\frac{d\bar{N}(t)}{dt} = \bar{A}\bar{N}(t)$$ where $\bar{A}$ is typically very sparse. The numerical difficulties associated with this set of coupled differential equations are twofold: