function [dist, generations] = simulate_population(starting_couples, chance_of_extinction)
% Sample input: [dist, generations] = simulate_population(2500, 1/2);
% evaluate dist(1) to get the probability of extinction after generations
transition_matrix = gen_transition_matrix(starting_couples);
dist = zeros(1, starting_couples + 1);
dist(end) = 1;
generations = 0;
while dist(1,1) <= chance_of_extinction
    dist = dist*transition_matrix;
    generations = generations + 1;
end

function transition_matrix = gen_transition_matrix(starting_couples)
transition_matrix = zeros(starting_couples + 1);
for i = 0:starting_couples
    binom_dist = makedist('Binomial', 'N', 2*i, 'p', 1/2);
    temp = pdf(binom_dist, 0:i);
    temp(1:end - 1) = 2*temp(1:end - 1);
    transition_matrix(i + 1, 1: i+1) = temp;
end