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

my_yprime = function(t, y, parameters){
    x = y[1]
    y = y[2]
    dy = numeric(2)
    dy[1] = 1
    dy[2] = ( (1+x)*y - 2*exp(-1/x) )/(x+1.e-3)^2 # the equation with tau=1/t
    list( dy )
}

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

diff_eq_params = c()

flowField(my_yprime, x.lim = c(1, 6), y.lim = c(-6, 4), parameters = diff_eq_params, points = 21, add = FALSE)
trajectory(my_yprime, y0 = c(1.0, -4.0), t.start=1.0, t.end = 6, parameters = diff_eq_params)
trajectory(my_yprime, y0 = c(1.0, -4/5), t.start=1.0, t.end = 6, parameters = diff_eq_params)
trajectory(my_yprime, y0 = c(1.0, exp(-1)), t.start=1.0, t.end = 6, parameters = diff_eq_params)
trajectory(my_yprime, y0 = c(1.0, 1.0), t.start=1.0, t.end = 6, parameters = diff_eq_params)

#dev.off()