# define proj matrix A
f2 <- 1
p1 <- 1

A <- rbind(c(0, f2),
           c(p1,0))

# Test for power-positivity.  n^2-2n+2 is crucial power.  For n=2, that's 2
A %*% A

eigen(A)

# Initial population size in each class
n_zero <- matrix((c(2900, 9000)), 2, 1)

Tmax <- 20
n_vs_t <- matrix(0, 2, Tmax)

n_vs_t[,1] <- n_zero

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

matplot(1:Tmax, t(n_vs_t), type = "l",
        xlab = "t", ylab = "n")
legend("topright", c("stage 1", "stage 2"), col = 1:2, lty = 1:2)