# iterates Leslie matrix

n_zero <- matrix(c(10, 100, 500), 3, 1)

A <- rbind(c(0, 1, 5),
           c(.5, 0, 0),
           c(0, .25, 0))

Tmax <- 200

# initialize a zero matrix of 3 rows and Tmax columns
n_vs_t <- matrix(0, 3, Tmax)

# put initial values of the population
n_vs_t[,1] <- n_zero

# %*% is matrix multiplication
for (t in 2:Tmax)
  n_vs_t[,t] <- A %*% n_vs_t[,t-1]

# plot two graphs on one graphics device
par(mfrow=c(2,1))
# matplot is for plotting matrices
matplot(1:Tmax, t(n_vs_t), type = "l", # t() is transpose of a matrix
        xlab = "t", ylab = "n", col = 1:3, lwd = 2, lty = 1) 
matplot(1:Tmax, log(t(n_vs_t)), type = "l", 
        xlab = "t", ylab = "log(n)", col = 1:3, lwd = 2, lty = 1) 

# fit a portion of log(n) to a straight line for all three populations
p <- matrix(0, 2, 3); a <- 25:Tmax
for (i in 1:3) {
  b <- log(n_vs_t)[i,25:Tmax]
  p[,i] <- lm(b ~ a)$coeff # $ extracts the coefficients of the linear fit
}

# a more parsimonious code that does the same thing
p <- apply(log(n_vs_t)[,25:Tmax], 1, function(x) lm(x ~ I(25:Tmax))$coeff)

lambda_estimate <- exp(p[2,1])

print('estimated lambda from matrix iteration is ...')
print(lambda_estimate)

#compute dominant eigenvalue of A
vA=eigen(A)

print('dominant eigenvalue is  ...')
 print(vA$values[1])