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

my_yprime = function(t, y, parameters){
    g = 9.8
    m = parameters[1]
    v_0 = parameters[2]
    gamma = g*m / v_0^2
    y_prime = g - ( gamma / m ) * y^2
    list( y_prime )
}

#postscript("../../WriteUp/Graphics/Chapter1/chap_1_sect_1_prob_25_plot.eps", onefile=FALSE, horizontal=FALSE)

diff_eq_params = c( 10, 49 ) # =( m, v_0 )

flowField(my_yprime, x.lim = c(0, 10), y.lim = c(40, 60), parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE)

nullclines(my_yprime, x.lim = c(0, 10), y.lim = c(40, 60), parameters = diff_eq_params, system ="one.dim")

trajectory(my_yprime, y0 = c(45.0, 55.0), t.end = 10, parameters = diff_eq_params, system = "one.dim")

#dev.off()