plot(1:days, R, pch=15, xlim=(c(0,days)), ylim=(c(0,ymax)), col="blue3", xlab="days", ylab="growth rate: N(t)/N(t-1)")

M <- matrix(0,nrow=stages,ncol=stages)		# The model is "V[,j+1] = M %*% V[,j]"

lines(1:days, R, pch=15, col="blue3", lty=2)

lambda1 <- Re(result$values[1])	# "Re(x)" ignores any imaginary part of "x" (should be 0)

s <- c(0.373,0.195,0.319,0.356,0.622,0.917)

N <- colSums(V)

for (t in 1:days)			# "t" is the new time value
	{	
	n <- t			# "n" is the column number for the old population values
	V[,n+1] <- M %*% V[,n]  	# "V[,n+1]" is the column of new population values
	}

V <- matrix(0,nrow=stages,ncol=days+1)		# "V[i,j]" is the population of stage "i" at time "j-1"

for (j in 1:(days+1)) 
	U[,j] <- V[,j]/N[j]

p <- c(0.593,0.805,0.681,0.489,0.358)

V0 <- c(0,0,0,0,0,1)

ymax <- min(max(R),3)

M[1,] <- f

abline(h=lambda1, lty=3)

result <- eigen(M)

V[,1] <- V0

R <- N[2:(days+1)]/N[1:days]	

proportions <- Re(result$vectors[,1]/sum(result$vectors[,1]))

N <- array(0,dim=days+1)			# "N[j]" is the total population at time "j-1"

M[6,4] <- a

R <- array(0,dim=days)			# "R[j]" is the ratio of "N" at time "j" to "N" at time "j-1"

for (j in 1:stages-1)
	M[j+1,j] <- p[j]

f <- c(0,0,0,0,0.9,4.2)

U <- matrix(0,nrow=stages,ncol=days+1)		# "U[i,j]" is the proportion of stage "i" at time "j-1"

r <- log(lambda1)			# "log" is natural logarithm

M[1,1] <- s[1]+f[1]
for (j in 2:stages)
	M[j,j] <- s[j]

days <- 20     

a <- 0.089

stages <- 6