# 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){
    x = y[1]
    y = y[2]
    dy = numeric(2)
    dy[1] = x - 4*y
    dy[2] = -x + y
    list(dy)
}

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

y_exact = function(t){
    x = (exp(-t) + exp(3*t))/2
    y = (exp(-t) - exp(3*t))/4
    c(x, y)
}

# Where do we want to evaluate the solution:
#
t_grid = seq( t0, t1, by=0.05)
t_grid = c( t1 ) 

# RK4 method for various h values: 
#
h = 0.2
n_iters = 6
sol_at_t1 = c()
for( ii in 1:n_iters ){
    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' )
    sol_at_t1 = rbind( sol_at_t1, c( h, rk_xhat$pred, rk_yhat$pred ) )

    h = h/2
}

R = data.frame(sol_at_t1)
colnames(R) = c('h', 'x_hat', 'y_hat')
print(R)

print(y_exact(1.0))