len <- function(v) sqrt(sum(v^2)) # Returns the length of v

horn <- matrix(c(0.05,0.36,0.5,0.09, 0.01,0.57,0.25,0.17,
            0,0.14,0.55,0.31, 0,0.01,0.03,0.96), 4,4) 
v0 <- c(1,0,0,0)
v1 <- horn %*% v0; v1
v2 <- horn %*% v1; v2

# Run n*50 years
v <- v0
data <- as.matrix(v)
n <- 15
for (i in seq(n)) {
  v <- horn %*% v
  data <- cbind(data,v)
}

# Now plot the data
time <- 50*seq(n,from=0)
plot(time,data[1,],xlab="Time in years",ylab="Frequency",type="b",pch=20,lwd=2,col="red")
lines(time,data[2,],type="b",pch=20,lwd=2,col="green")
lines(time,data[3,],type="b",pch=20,lwd=2,col="blue")
lines(time,data[4,],type="b",pch=20,lwd=2,col="gold")
legend("topright",legend=c("Gray birch","Blackgum","Red maple","Beech"),lty=1,col=c("red","green","blue","gold"),cex=0.75)

# The climax is given by the dominant eigenvector:
e <- eigen(horn)
ev1 <- e$vectors[,1]; ev1
ev1 <- ev1/sum(ev1); ev1
data[,n]

# The proper way to do this in R is not via a loop but by using apply:
n <- 25
v <- v0
data <- sapply(seq(n),function(i){v <<- horn %*% v})
data <- cbind(as.matrix(v0),data)
data[,n]

# Make figure
size <- 5 
pdf("horn.pdf",width=size,height=size)
#opar <- par(no.readonly=TRUE) 
par(mar=c(2.6,2.6,1.2,0.2),mgp=c(1.5,0.5,0)) 
time <- 50*seq(n,from=0)
plot(time,data[1,],xlab="Time in years",ylab="Frequency",type="b",pch=20,lwd=2,col="red")
lines(time,data[2,],type="b",pch=20,lwd=2,col="darkgreen")
lines(time,data[3,],type="b",pch=20,lwd=2,col="blue")
lines(time,data[4,],type="b",pch=20,lwd=2,col="gold")
legend("topright",legend=c("Gray birch","Blackgum","Red maple","Beech"),lty=1,col=c("red","darkgreen","blue","gold"),lwd=2,cex=0.75)
mtext("(a)",3,0.1,at=350)
dev.off()