if( !require('phaseR') ){
    install.packages('phaseR')
}
library(phaseR)

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    dy = rep(NA, length(y))
    dy[1] = xt*(1.0 - parameters$sigma * xt - 0.5 * yt)
    dy[2] = yt*(-0.75 + 0.25 * xt)
    list(dy)
}

t.end = 30.0

#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_5_prob_11_plot.eps', onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1, 2))

# One value of sigma:
# 
diff_eq_params = list(sigma=0.1)
CP = data.frame(
    x=c(0, 1/diff_eq_params$sigma, 3),
    y=c(0, 0, 2-6*diff_eq_params$sigma)
)
print(CP)

flowField(my_yprime, x.lim = c(-0.1, 12.0), y.lim = c(-0.1, 4.5), parameters = diff_eq_params, points = 30, add = FALSE, xlab='x', ylab='y', main=sprintf('sigma= %f', diff_eq_params$sigma))
trajectory(my_yprime, y0 = c(0.1, 4), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)


# A second value of sigma:
# 
diff_eq_params = list(sigma=(0.2637626 + 0.01))
CP = data.frame(
    x=c(0, 1/diff_eq_params$sigma, 3),
    y=c(0, 0, 2-6*diff_eq_params$sigma)
)
print(CP)

flowField(my_yprime, x.lim = c(-0.1, 12.0), y.lim = c(-0.1, 4.5), parameters = diff_eq_params, points = 30, add = FALSE, xlab='x', ylab='y', main=sprintf('sigma= %f', diff_eq_params$sigma))
trajectory(my_yprime, y0 = c(0.1, 4), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

par(mfrow=c(1, 1))

#dev.off()