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

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    zt = y[3]
    dy = rep(NA, length(y))
    dy[1] = parameters$sigma * (-xt + yt)
    dy[2] = parameters$r * xt - yt - xt*zt
    dy[3] = -parameters$b * zt + xt*yt
    list(dy)
}

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

diff_eq_params = list(sigma=10, b=8/3, r=28)

ts = seq(0, 20, length.out=2000)

res_1 = ode(y=c(5, 5, 5), times=ts, fun=my_yprime, parms=diff_eq_params)
xs = res_1[, 2] # the numerical function x(t)
ys = res_1[, 3] # the numerical function y(t)
zs = res_1[, 4] # the numerical function z(t)

par(mfrow=c(1, 2))

plot(xs, ys, type='l', xlab='x(t)', ylab='y(t)')
grid()

plot(xs, zs, type='l', xlab='x(t)', ylab='z(t)')
grid()

par(mfrow=c(1, 1))

#dev.off()