source('../Chapter2/euler.R') source('improved_euler.R') t0 = 0.0 y0 = 1 dy.dt = function(t, y){ 0.5 - t + 2*y } h=0.1 e_res = euler(dy.dt, h=h, start=t0, y0=y0, end=t0+h) ie_res = improved_euler(dy.dt, h=h, start=t0, y0=y0, end=t0+h) # The estimated error: # e_est = abs(ie_res$ys[2] - e_res$ys[2]) e_target = 0.0025 factor = sqrt(e_target/e_est) print(sprintf('h= %.2f; ieuler est= %f; euler est= %f; error est= %f; factor= %f; new_h= %f', h, ie_res$ys[2], e_res$ys[2], e_est, factor, factor*h))