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

A = matrix(c(1/4, -2, 1, -5/4), nrow=2, ncol=2, byrow=T)
print(eigen(A))

# Part (b) Plot (x1, x2):
#
my_yprime = function(t, y, parameters){
    dy = parameters$A %*% y
    list(dy)
}
diff_eq_params = list(A=A)

#postscript("../../WriteUp/Graphics/Chapter7/chap_7_sect_6_prob_11_x1_x2_plot.eps", onefile=FALSE, horizontal=FALSE)
flowField(my_yprime, x.lim=c(-10, 10), y.lim=c(-10, 10), parameters = diff_eq_params, points=21, add=FALSE, xlab='x_1', ylab='x_2', main='Problem 11')
trajectory(my_yprime, y0 = c(3, 8), t.end = 3.0, parameters = diff_eq_params, col='black')
trajectory(my_yprime, y0 = c(-2.5, -5), t.end = 2.0, parameters = diff_eq_params, col='black')
trajectory(my_yprime, y0 = c(5, -5), t.end = 1.0, parameters = diff_eq_params, col='black')
trajectory(my_yprime, y0 = c(-2, 2), t.end = 1.0, parameters = diff_eq_params, col='black')
grid()
#dev.off()


# Part (b) Plot (t, x1(t)) and (t, x2(t)):
#
t = seq( 0, 15, length.out=1000 )

state = c( y=c(3, 8) )

parms = list(A=A)
res = ode(y=state, times=t, fun=my_yprime, parms=parms)
ts = res[,'time']

#postscript("../../WriteUp/Graphics/Chapter7/chap_7_sect_6_prob_11_t_x1_N_t_x2.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1, 2))
plot(ts, res[,'y1'], type='l', lty=1, xlab='t', ylab='y_1(t)')
lines(c(0), c(3), type='p', pch=19, cex=1.5)
grid()
plot(ts, res[,'y2'], type='l', lty=1, xlab='t', ylab='y_2(t)')
lines(c(0), c(8), type='p', pch=19, cex=1.5)
grid()
par(mfrow=c(1, 1))
#dev.off()