% program to simulate a first-order reaction
clear % clear the stored variables
Nmol=100; % number of molecules
x=ones(Nmol,1); % initial state of molecules is 1
k=0.05;% first-order rate constant (s-1)
tau=1000;% time step (s)
Nsteps=100; % number of time steps in the simulation
p=1-exp(-k*tau); % probability of the reaction occurring in time tau
t(1)=0; % initial time (s)
y(1)=Nmol; % initial number of molecules in state 1
for j=1:Nsteps-1 % take Nsteps-1 steps
    for i=1:Nmol % go through each molecule and see if it reacted
        tmp=rand; % generate a random number
        if tmp<p  % compare the random number to prob of rxn
            x(i)=0; % if rand#<prob, then change state to 0
        end
    end
    yytmp=mean(x); % calculate mean state of molecules
    y(j+1)=Nmol*yytmp; % record the number of molec in state 1 
    t(j+1)=t(j)+tau; % increment the time
end
plot(t,y) % plot the mean state as a function of time
axis([0  100 0 Nmol]) % adjust the axes on the plot
xlabel('TIME (s)','FontSize',16)
ylabel('NUMBER OF MOLECULES','FontSize',16)