function [t,n] = birthdeath (beta, mu, K, n0, T)

% BIRTHDEATH Nonlinear stochastic birth death process simulation.
%             [t,n] = birthdeath (beta, mu, K, n0, T)
%            Input:
%             beta = per capita birth rate
%             mu   = per capita baseline death rate
%             K    = "carrying capacity"
%             n0   = initial population size
%             T    = number of steps to take (default = 100)
%            Output:
%             t    = times at which events (births or deaths) occur
%             n    = corresponding population size
%           

if (nargin < 5), T = 100; end           % if T is not given, set T = 100

n = zeros(T+1,1);
t = zeros(T+1,1);

n(1) = n0;
t(1) = 0;

for i = 1:T,
  if (n(i) == 0),                       % if n = 0, truncate t & n, and loop no more
    n = n(1:i);
    t = t(1:i);
    break;
  end
  brate = beta * n(i);                  % population birth rate
                                        % population death rate
  drate = (mu + (beta - mu)*n(i)/K) * n(i);
  erate = brate + drate;                % total event rate
  wt = -log(rand)/(erate);              % exponentially-distributed waiting time
  dn = 2*(rand < brate/erate)-1;        % is it a birth or death?
  t(i+1) = t(i) + wt;
  n(i+1) = n(i) + dn;
end