# Draw the phase plane of a linear 2D system.

# The system
A <- 0.5 * rbind( c(1, -3), c(-3, 1) )

# Eigenvectors
ev1x <- c(-4,4)
ev1y <- c(-4,4)
ev2x <- c(-4,4)
ev2y <- c(4,-4)
matplot( ev1x, ev1y, xlab = "x1", ylab = "x2", type = "l",col="red", xlim=c(-4,4), ylim=c(-4,4), add = FALSE )
matplot( ev2x, ev2y, type = "l",col="red", add = TRUE )

# Nullclines
nc1x <- c(-6,6)
nc1y <- c(-2,2)
nc2x <- c(-2,2)
nc2y <- c(-6,6)
matplot( nc1x, nc1y, type = "l",col="green", add = TRUE )
matplot( nc2x, nc2y, type = "l",col="green", add = TRUE )

# Trajectories
maxTime <- 10
dt <- 0.01 # Precision of simulation
N <- maxTime / dt # Number of samples
X <- matrix( 0, nrow=2, ncol=N ) # The results will go here

# Run the simulation
for( j in 1:6 ) {
  xy <- locator( n=1 ) # get initial position from mouse click
  X[1,1] <- xy[[1]]
  X[2,1] <- xy[[2]]
  for( t in 2:N ) {
    dx <- A %*% X[1:2,t-1]
    X[1:2,t] <- X[1:2,t-1] + dt*dx
  }
  matplot(X[1,1:N],X[2,1:N],type = "l",col="black", add = TRUE )
}