# 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, u, parameters){
    with(as.list(c(u, parameters)),{
        t = u[1]
        y = u[2]
        du = numeric(2)
        du[1] = 1
        du[2] = lambda * y + 1 - lambda * t
        list(du)
    })
}

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

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

# Where will we compute the solution:
#
h = 0.01

compute_grid = seq( t0, t1, by=h )

# Integrate the ODE for each value of lambda:
#
lambdas = c(1, 10, 20, 50)
YHAT = c()
for( lam in lambdas ){
    res = ode(y0, compute_grid, my_yprime, parms=c(lambda=lam), method='rk4')
    res = as.data.frame(res)
    colnames(res) = c('time','x','y')
    yhat = knn.reg( as.matrix(res$time), y=res$y, test=as.matrix(t_grid), k=1, algorithm='brute')
    YHAT = cbind( YHAT, yhat$pred )
}
YHAT = cbind(t_grid, YHAT)
colnames(YHAT) = c('t', sprintf('lam= %.0f', lambdas) )
print(YHAT)