"birthdeath"<-
function(a = 0.1, lambda = 0.1, mu = 0.05, N = 1000, initial.size = 10)
{
	#
	# To set the random number seed, un-comment the following line and 
	# plug in any non-negative integer in [0,1023].
	# set.seed(123)
	#
	# Simulate N events from a linear birth and death process with immigration.
	# a = immigration rate
	# 
	S <- rep(0, N)
	W <- S
	X <- S
	X[1] <- initial.size
	W[1] <- S[1]
	# Note that this is time = 0.
	#
	for(i in 2:N) {
		#
		# Simulate a sojourn (inter-event) time, and compute the waiting time.
		#
		if(a>0 | X[i-1]>0){
			S[i] <- rexp(1, rate = a + (lambda + mu) * X[i - 1])
			W[i] <- W[i - 1] + S[i]
			#
			# Determine the probability that the event is a death.
			#
			death <- (mu * X[i - 1])/(mu * X[i - 1] + a + lambda * X[i -
				1])
			#
			# Sample from the set of possible outcomes: -1 for death or +1 for birth. 
			#
			X[i] <- X[i - 1] + sample(c(-1, 1), size = 1, prob = c(death,
				1 - death))
			}#end if
		else{
			if(a>0){#immigration saves the extinct population
				S[i]_rexp(1,rate=a)
				W[i]_W[i-1]+S[i]
				X[i]_X[i-1]+1
				}# end if a>0
			else{#population is extinct, with no hope of immigration
				S[i]_0
				W[i]_W[i-1]
				X[i]_X[i-1]
				}
		}#end else
	}
	# end for loop on i
	#
	# Compute the analytic mean function of the process on a fine grid of times.
	# 
	grid <- ((0:1000) * max(W))/1000
	M <- (a/(lambda - mu)) * (exp((lambda - mu) * grid) - 1) + initial.size *
		exp((lambda - mu) * grid)
	#
	# Set up a plotting region, then add curve for the analytical mean and step function
	# for the birth and death process.
	#
	plot(c(W, grid), c(X, M), type = "n", xlab = "Time", ylab = 
		"Population Size")
	lines(grid, M, col = 8)
	lines(W,X,type="s")
	return()
}