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_6_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_1 = res_1[, 2] # the numerical function x(t)

res_2 = ode(y=c(5.01, 5, 5), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_2 = res_2[, 2] # the numerical function x(t)

Xs = cbind(xs_1, xs_2)

matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)')
grid()

#dev.off()