function [Svec, Ivec] = simDisease(contacts_per_day, prob_recovery, r_frac)
% Run a discrete SIR model to simulate the spread of a disease.
% `contacts_per_day` is the number of contacts people have each day.
% `prob_recovery` is the probability that an infectious person recovers
% each day. `r_frac` (betwteen 0 and 1) is the initial fraction of the
% population that is already resistant to the disease.
% Returns two vectors containing the number of susceptible (`Svec`) and
% infectious (`Ivec`) people on each day. Index 1 corresponds to day 0
% (initial conditions).
% The total population size is 10000, the initial number of infected people
% is 100, and the probability of infection during a contact is 30%. The
% model is run for 30 days (so the output vectors will contain 31
% elements).