# 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)

# Use the ode function in the deSolve library to integrate the differential equations:
# 
if( !require('deSolve') ){
    install.packages('deSolve')
}
library(deSolve)

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

t0 = 0.0
t1 = 1.0
y0 = c(1.0, 2.0)

# Where do we want to evaluate the solution:
#
t_grid = c( 0.5, 1.0 ) 

# RK4 method for various h values: 
#
h = 0.1
compute_grid = seq( t0, t1, by=h )
res = ode(y0, compute_grid, my_yprime, parms=c(), method='rk4')
res = as.data.frame(res)
colnames(res) = c('time', 'x', 'y')

rk_xhat = knn.reg( as.matrix(res$time), y=res$x, test=as.matrix(t_grid), k=1, algorithm='brute' )
rk_yhat = knn.reg( as.matrix(res$time), y=res$y, test=as.matrix(t_grid), k=1, algorithm='brute' )

R = data.frame(rk_xhat$pred, rk_yhat$pred)
R = cbind(t_grid, R)
colnames(R) = c('time', 'x', 'xprime')
print(R)