% Setting Initial Parameters
N = 2000;
init_pop = [N, 0, 0, 0];
half_lives = [2.5e5, 8e4, 1.6e3, 4/365];
t = linspace(0, 1e6, 1000);

%Calling monte carlo simulation function
pop_time_data = simulate_decay ( init_pop , half_lives , t);

l = log(2)./half_lives;
% Determining analytical solution
N1 = N*exp(-l(1)*t);
N2 = N*(l(1)/(l(2)-l(1)))*(exp(-l(1)*t)-exp(-l(2)*t));
N3 = N*l(1)*l(2)*((exp(-l(1)*t)/((l(1)-l(2))*(l(1)-l(3))))+(exp(-l(2)*t)/((l(2)-l(3))*(l(2)-l(1))))+(exp(-l(3)*t)/((l(3)-l(1))*(l(3)-l(2)))));
N4 = N - (N1+N2+N3);

hold on

%Plotting analytical solution
plot(t,N4,'k')
plot(t,N3,'g')
plot(t,N2,'r')
plot(t,N1,'b')

%Plotting monte carlo solution data
plot(t, pop_time_data(:,1)','--')
plot(t, pop_time_data(:,2)','--')
plot(t, pop_time_data(:,3)','--')
plot(t, pop_time_data(:,4)','--')

xlabel("Time (years)")
ylabel("Number of particles")
hold off