% EULERTAYLOR small code used to generate numbers in Dec 14 class

format short;
f = @(t,x) 1.025 * x;
t0 = 0;
x0 = 1;
tf = 3;

exact = exp(1.025 * tf);

for n=10.^(1:5)
  t = t0;  x = x0;                      % reset to initial condition
  h = (tf - t)/n;
  % Euler's method
  for j=1:n
    x = x + h * f(t,x);
    t = t + h;
  end
  printf('Euler error for     n = %7d:  %.3e\n',n,abs(x-exact))
end

f_xp = @(t,x) 1.025 * x;                % same as f, re-written with notation xp = x'
df_xpp = @(t,x,xp) 1.025 * xp;          % computes xpp = x''
ddf_xppp = @(t,x,xp,xpp) 1.025 * xpp;    % computes xppp = x'''
for n=10.^(1:3)
  t = t0;  x = x0;                      % reset to initial condition
  h = (tf - t)/n;
  % Taylor order 3
  for j=1:n
    xp = f_xp(t,x);
    xpp = df_xpp(t,x,xp);
    xppp = ddf_xppp(t,x,xp,xpp);
    x = x + h * (xp + (h/2) * (xpp + (h/3) * xppp));  % use Horner to eval poly in h (!)
    t = t + h;
  end
  printf('Taylor O3 error for n = %7d:  %.3e\n',n,abs(x-exact))
end
