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] = -yt - zt
    dy[2] = xt + parameters$a * yt
    dy[3] = parameters$b + zt * ( xt - parameters$c )
    list(dy)
}

my_jacobian = function(x, y, z, c){
    ## returns the Jacobian at the point(x, y, z)
    ##
    J = matrix(c(0, -1, -1, 1, 1/4, 0, z, 0, x-c), nrow=3, ncol= 3, byrow=TRUE)
}

# Parts (a):
#

## c=3:
##
diff_eg_params = list(a=1/4, b=1/2, c=3)
zs = 2*diff_eq_params$c + 2*sqrt(diff_eq_params$c*2 - 0.5) * c(-1, +1)
xs = diff_eq_params$a * zs
ys = -zs
CP = data.frame(x=xs, y=ys, z=zs)

As = mapply(my_jacobian, CP$x, CP$y, CP$z, diff_eq_params$c, SIMPLIFY=FALSE)
for(ii in 1:length(As) ){
    es = eigen(As[[ii]])
    print(sprintf('CP: x= %f; y= %f; z= %f', CP$x[ii], CP$y[ii], CP$z[ii]))
    print('Jacobian=')
    print(As[[ii]])
    print('eigenvalues=')
    print(round(es$values, 4))
}

## Part (b):
##
ts = seq(0, 20, length.out=2000)

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

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


## Part(c):
##
## Plot the solutions:
##
par(mfrow=c(3, 1))

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

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

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

par(mfrow=c(1, 1))