number_simulations := 1000:
initial_population_size := 100:
probability_bacterium_dies := .5:
T := 1000: # maximum number of generations allowed
generation_of_extinction := Array(1..number_simulations): # initialized to zero
population_size := Array(1..number_simulations, 1..T);
for i from 1 to number_simulations by 1 do
generation := 0:
pop_size := initial_population_size:
while (pop_size > 0 and generation < T) do
generation := generation + 1; 
x := pop_size; # we need the loop below to have a fixed limit - pop_size will be varying
for j from 1 to x by 1 do # determine for each bacterium whether it dies or reproduces
u := rand()/10^12;
if (u <= probability_bacterium_dies) then pop_size := pop_size-1 else pop_size := pop_size + 1 end if;
end do:
population_size[i, generation] := pop_size;
end do:
          # record the generation of extinction or T if it is smaller
generation_of_extinction[i] := generation; 
end do:

#############################################################
## The following code is how I wrote the data in the 
## generation_of_extinction array to a file. You will have to
## change the path, the 
## C://Users//Elizabeth Housworth//Desktop//generations_1.txt
## part. 

writedata("C://Users//Elizabeth Housworth//Desktop//generations_1.txt", convert(generation_of_extinction, list));
fclose("C://Users//Elizabeth Housworth//Desktop//generations_1.txt");

## The following code is how I wrote the data in the 
## population_size array to a file. You will have to
## change the path, the 
## C://Users//Elizabeth Housworth//Desktop//popsize_1.txt
## part. 

writedata("C://Users//Elizabeth Housworth//Desktop//popsize_1.txt", convert(population_size, listlist));
fclose("C://Users//Elizabeth Housworth//Desktop//popsize_1.txt");

###########################################################
## After you have written data to a file, you might want to get it 
## back sometime. To do that, you can use the following command
## which will open a file dialog box from which you can choose
## the file containing your data. 

ImportData();



#############################################################
## The following code uses the statistics package to calculate
## the mean population size and its standard error for each generation


with(Statistics):
average_number_bacteria := Array(1..1000):
SE_bacteria := Array(1..1000):
for l from 1 to 1000 by 1 do
average_number_bacteria[l] := Mean(population_size[1..1000,l]);
SE_bacteria[l] := StandardError(Mean, population_size[1..1000,l]);
end do:


##############################################################
## The following code creates an error plot (means with standard
## error bars) for the average number of bacteria in each generation
## It uses only every 10th generation - otherwise the plot is too
## cluttered

with(Statistics):
ErrorPlot(average_number_bacteria[[seq(10*i, i=1..100)]], yerrors=SE_bacteria[[seq(10*i, i=1..100)]]);



####################################################
# vary probability bacterium dies each generation

number_simulations := 1000:
initial_population_size := 100:
T := 1000: # maximum number of generations allowed
generation_of_extinction := Array(1..number_simulations): # initialized to zero
population_size := Array(1..number_simulations, 1..T);
for i from 1 to number_simulations by 1 do
generation := 0:
pop_size := initial_population_size:
while (pop_size > 0 and generation < T) do
generation := generation + 1; 
probability_bacterium_dies := .4 + evalf(rand()/(5*10^12)):
x := pop_size; # we need the loop below to have a fixed limit - pop_size will be varying
for j from 1 to x by 1 do # determine for each bacterium whether it dies or reproduces
u := rand()/10^12;
if (u <= probability_bacterium_dies) then pop_size := pop_size-1 else pop_size := pop_size + 1 end if;
end do:
population_size[i, generation] := pop_size;
end do:
          # record the generation of extinction or T if it is smaller
generation_of_extinction[i] := generation; 
end do:


#############################################################
## The following code is how I wrote the data in the 
## generation_of_extinction array to a file. You will have to
## change the path, the 
## C://Users//Elizabeth Housworth//Desktop//generations_2.txt
## part. 

writedata("C://Users//Elizabeth Housworth//Desktop//generations_2.txt", convert(generation_of_extinction, list));
fclose("C://Users//Elizabeth Housworth//Desktop//generations_2.txt");

## The following code is how I wrote the data in the 
## population_size array to a file. You will have to
## change the path, the 
## C://Users//Elizabeth Housworth//Desktop//popsize_2.txt
## part. 

writedata("C://Users//Elizabeth Housworth//Desktop//popsize_2.txt", convert(population_size, listlist));
fclose("C://Users//Elizabeth Housworth//Desktop//popsize_2.txt");