source('./improved_euler.R')

t0 = 0.0
y0 = 1.0 
t1 = 4
h = 0.05

# Where do we want to evaluate the solution:
#

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

# The exact solution:
#
y_exact_fn = function(t){ (4*t-3)/16 + (19/16)*exp(4*t) }
y_exact = y_exact_fn(x_grid)

# Compute the local truncation error bound for the improved Euler method: 38/3) e^{4 \bar{t}_n} h^3:
#
lteb = (38/3) * exp( 4*(x_grid + h) ) * h^3 # holds |e_1|, |e_2|, |e_3|, |e_4|
lteb = lteb[c(1, 4)]

# Compute the two true errors at y_1 and y_4:
#
dy.dt = function(t, y){ 1 - t + 4*y }

res = improved_euler(dy.dt, h=h, start=t0, y0=y0, end=t0+h)
e_y1 = res$ys[2]

t3 = t0 + 3*h
y3 = y_exact_fn(t3)
res = improved_euler(dy.dt, h=h, start=t3, y0=y3, end=t3+h)
e_y4 = res$ys[2]

# The true error for that improved Euler step:
#
true_error = y_exact[c(1, 4)+1] - c( e_y1, e_y4 )
A = rbind( lteb, abs(true_error) )
rownames(A) = c('Local Truncation Error:', 'abs(True Error):')
colnames(A) = c( 'abs(e1)', 'abs(e4)' )
print(A)