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

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

my_jacobian = function(x, y) {
    # returns the Jacobian at the point (x, y)
    #
    J = matrix(c(-2 - 3*x^2 - y^2, -1 - 2*x*y, 1 + 2*x*y, -1 + x^2 + 3*y^2), nrow=2, ncol=2, byrow=TRUE)
}

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    dy = rep(NA, length(y))
    dy[1] = -2*xt - yt - xt*(xt^2 + yt^2)
    dy[2] = xt - yt + yt*(xt^2 + yt^2)
    list(dy)
}

diff_eq_params = list()

print('Solving f(x) = 0 numerically to find the roots:')

rts_1 = nleqslv(c(1, 1), function(x){my_yprime(NA, x, diff_eq_params)[[1]]}, jac=function(x){my_jacobian(x[1], x[2])})
print(rts_1$x)

rts_2 = nleqslv(c(-1, -1), function (x){my_yprime(NA, x, diff_eq_params)[[1]]}, jac=function(x){my_jacobian(x[1], x[2])})
print(rts_2$x)

# Classify the critical points of this system:
CP = data.frame (x=c(rts_1$x[1], 0, rts_2$x[1]), y=c(rts_1$x[2], 0, rts_2$x[2]))

As = mapply(my_jacobian, CP$x, CP$y, SIMPLIFY=FALSE)

for (ii in 1:length(As)){
    es = eigen(As[[ii]])
    print(sprintf('CP: x= %f; y= %f', CP$x[ii], CP$y[ii]))
    print('Jacobian=')
    print(As[[ii]])
    print('eigenvalues=')
    print(es$values)
}

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

t.end = 5.0
L = 4
flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y')

# Plot some trajectories with initial condition near the critical points:
#
trajectory(my_yprime, y0 = c(1, 1), t.end = t.end, parameters = diff_eq_params, col='black', pch=19)
trajectory(my_yprime, y0 = c(-1, 1), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
trajectory(my_yprime, y0 = c(1, -1), t.end = t.end, parameters = diff_eq_params, col='red', pch=19)
trajectory(my_yprime, y0 = c(-1, -1), t.end = t.end, parameters = diff_eq_params, col='green', pch=19)

#dev.off()