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

m <- matrix(c(1,2,2,1),2,2)

v0 <- c(0,1)
v1 <- m %*% v0; v1
v2 <- m %*% v1; v2
v3 <- m %*% v2; v3

v <- v0
plot(c(0,v[1]),c(0,v[2]),xlab="x",ylab="y",type="l",xlim=c(0,15),ylim=c(0,15),lwd=2,col=2)
points(v[1],v[2],pch=20,col=2)
text(v[1]+0.5,v[2],"0")
for (i in seq(3)) {
  v <- m %*% v; print(c(i,v))
  lines(c(0,v[1]),c(0,v[2]),lwd=2,col=2+i)
  points(v[1],v[2],col=2+i,pch=20)
  text(v[1]+0.5,v[2],i)
}

v <- v0 
s <- v/len(v)
plot(c(0,s[1]),c(0,s[2]),xlab="x",ylab="y",type="l",xlim=c(0,1),ylim=c(0,1),col=2)
points(s[1],s[2],col=2)
text(s[1]+0.05,s[2],"0")
for (i in seq(5)) {
  v <- m %*% v; print(c(i,v)); s <- v/len(v)
  lines(c(0,s[1]),c(0,s[2]),col=2+i)
  points(s[1],s[2],col=2+i)
  text(s[1]+0.05,s[2],i)
}

eigen(m)

# Make figure
size <- 5 
pdf("converge.pdf",width=size,height=size)
par(mar=c(2.6,2.6,1.2,0.2),mgp=c(1.5,0.5,0)) 
v <- v0
plot(c(0,v[1]),c(0,v[2]),xlab="x",ylab="y",type="l",xlim=c(0,15),ylim=c(0,15),lwd=2,col="red")
points(v[1],v[2],pch=20,col="red")
v <- m %*% v
lines(c(0,v[1]),c(0,v[2]),lwd=2,col="darkgreen")
points(v[1],v[2],pch=20,col="darkgreen")
v <- m %*% v
lines(c(0,v[1]),c(0,v[2]),lwd=2,col="blue")
points(v[1],v[2],pch=20,col="blue")
v <- m %*% v
lines(c(0,v[1]),c(0,v[2]),lwd=2,col="gold")
points(v[1],v[2],pch=20,col="gold")
legend("topleft",legend=c("(0,1)","(2,1)","(4,5)","(14,13)"),lty=1,col=c("red","darkgreen","blue","gold"),lwd=2)
mtext("(b)",3,0.1,at=7)
dev.off()