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

if( !require('FNN') ){
    install.packages("FNN", dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(FNN)

source('euler.R')

my_yprime = function(t, y, parameters){
    x = y[1]
    y = y[2]
    dy = numeric(2)
    dy[1] = 1
    dy[2] = -x*y + 0.1*y^3
    list( dy )
}

#postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_7_prob_18_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c()

logistic.flowField <- flowField(my_yprime, x.lim=c(0, 5), y.lim=c(-5, 5), parameters=diff_eq_params, points=21, add=FALSE, main='Problem 18')

logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 2.5 ), t.end = 3.0, parameters = diff_eq_params, col='black') # diverges

logistic.trajectory <- trajectory(my_yprime, y0 = c( 0, 2.0 ), t.end = 3.0, parameters = diff_eq_params, col='black') # converges

#dev.off()

# Thus we know that 2.0 (converges) < alpha_0 < 2.5 (diverges)
# Use a bisection like algorithm to reduce the size of the interval.
#
bkt = c( 2.0, 2.5 )
for( iter in 1:10 ){
    alpha_test = mean(bkt)
    dy.dx = function(x, y){ -x*y + 0.1*y^3 }
    h = 0.01
    res = euler(dy.dx, h=h, start=0, y0=alpha_test, end=3.0)
    ys = res$ys
    ys = ys[is.finite(ys)]
    npts = length(ys)
    if( ys[npts] < 0.5 ){ # converges
        bkt[1] = alpha_test
    }else{ # diverges
        bkt[2] = alpha_test
    }
    print( bkt )
}