model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dx <- a*x + b*y
    dy <- c*x + d*y
    return(list(c(dx, dy)))  
  }) 
}  

algebra <- function(parms) {
  with(as.list(parms), {
    tr <- a+d
    det <- a*d - b*c
    D <- tr^2 - 4*det
    cat("Trace: ",tr,"Determinant:",det,"Discrimant:",D,"\n")
    if (D>0) {
      l1 <- (tr+sqrt(D))/2
      l2 <- (tr-sqrt(D))/2
      cat("Eigenvalues:",l1,l2,"Eigenvectors: (",-b,a-l1,"), (",-b,a-l2,")")
    } else {
      re <- tr/2
      im <- sqrt(-D)/2
      cat("Complex eigenvalues:",re,"+/-",im,"i")
    }
  })
}

s <- c(x=0,y=0)

# To make a stable node we need tr<0, det>0, and D>0
# we therefore make the diagonal elements larger than the off-diagonal elements:
p <- c(a=-2,b=1,c=1,d=-2)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# To make an unstable node we need tr>0, det>0, and D>0
# so we reverse the sign of the large diagonal elements:
p <- c(a=2,b=1,c=1,d=2)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# To make a saddle point we need det<0, so we need small diagonal elements
p <- c(a=-1,b=-2,c=-2,d=-1)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# To make a stable spiral point we need tr<0, det>0, and D<0 
# so we need small diagonal elements and make a large determinant by setting c=-2
p <- c(a=-1,b=2,c=-2,d=-1)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# To make an unstable spiral point we need tr>0, det>0, and D<0 
# so we reverse the sign of diagonal elements
p <- c(a=1,b=2,c=-2,d=1)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# To make a center point we need tr=0, det>0, and D<0 
# so we need small diagonal elements and make a large determinant by setting c=-2
p <- c(a=0,b=2,c=-2,d=0)
algebra(p)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10)
newton(s,jacobian=TRUE,vector=TRUE)

# Make figure
size <- 4 
pdf("linear.pdf",width=2*size,height=3*size)
opar <- par(no.readonly=TRUE) 
par(mar=c(2.6,2.6,1.2,0.2),mgp=c(1.5,0.5,0)) 
par(mfrow=c(3,2)) 
p <- c(a=-2,b=1,c=1,d=-2)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(a)")
p <- c(a=2,b=1,c=1,d=2)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(b)",legend=FALSE)
p <- c(a=-1,b=-2,c=-2,d=-1)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(c)",legend=FALSE)
p <- c(a=-1,b=2,c=-2,d=-1)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(d)",legend=FALSE)
p <- c(a=1,b=2,c=-2,d=1)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(e)",legend=FALSE)
p <- c(a=0,b=2,c=-2,d=0)
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=T,tstep=0.1,tmax=10,main="(f)",legend=FALSE)
dev.off()

size <- 5 
pdf("stableNode.pdf",width=size,height=size)
opar <- par(no.readonly=TRUE) 
par(mar=c(3.6,2.6,1.2,0.2),mgp=c(1.5,0.5,0)) 
p <- c(a=-2,b=1,c=1,d=-2)
colors[1] <- "white"; colors[2] <- "white"
plane(xmin=-1,xmax=1,ymin=-1,ymax=1,portrait=TRUE,tstep=0.1,tmax=10,lwd=0,legend=FALSE)
lines(c(-1,1),c(-1,1),col="red",lwd=2)
lines(c(-1,1),c(1,-1),col="blue",lwd=2)
lines(c(0,0),c(-1,1))
lines(c(-1,1),c(0,0))
run(10,0.1,state=c(x=1,y=0.5),col="darkgreen",traject=TRUE)
dev.off()