# Part (a):
#
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] = x^2 + y^2
    list(dy)
}

#postscript("../../WriteUp/Graphics/Chapter8/chap_8_sect_3_prob_14_direction_field.eps", onefile=FALSE, horizontal=FALSE)
diff_eq_params = c()
flowField(my_yprime, x.lim = c(0, 1.1), y.lim=c(-1, 4), parameters = diff_eq_params, points = 21, add = FALSE, xlab = 't')
trajectory(my_yprime, y0 = c(0.0, 1.0), t.end=0.99, parameters= diff_eq_params)
#dev.off()


# Part (b):
#
# I use the function knn.reg to evaluate the numerical solution at a fixed grid of points
#
if( !require('FNN') ){
    install.packages('FNN', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(FNN)

source('runge_kutta.R')

t0 = 0.0
y0 = 1
t1 = 1 #0.95

# Where do we want to evaluate the solution:
#
x_grid = c( 0.8, 0.9, 0.95, 1.0 )

dy.dt = function(t, y){ t^2 + y^2 }

h=0.01
res = runge_kutta(dy.dt, h=h, start=t0, y0=y0, end=t1)
rk_yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

A = rbind( rk_yhat_h1$pred )
rownames(A) = NULL
R = cbind( c( 0.05 ), A )
colnames(R) = c( 'h', sprintf( 'y(%.2f)', x_grid ) )
print(R)