# 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('improved_euler.R')

t0 = 0.0
y0 = 1.0
t1 = 2.0

# Where do we want to evaluate the solution:
#
x_grid = c( seq( t0, 0.5, by=0.1 ), 1.0, 1.5, 2.0 ) # TODO: fix limits

# The improved Euler method:
#
dy.dt = function(t, y){ 1-t+4*y }

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

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

A = rbind( ie_yhat_h1$pred, ie_yhat_h2$pred )
rownames(A) = NULL
R = cbind( c( 0.025, 0.01 ), A ) # TODO: should be 0.1
colnames(R) = c( 'h', sprintf( 'y(%.2f)', x_grid ) )
print('Problem 13 (Improved Euler):')
print(R)