% Solves IVP in 5.2.1b
f = @(t,y) 1 + (t-y)^2;
df = @(t,y) 2*(t-y) - 2*f(t,y)*(t-y);
y = @(t) t + 1/(1-t); % true solution
a=2; b=3; n=2; y0=1; h=(b-a)/n;

y1 = euler(a,b,n,y0,f);
y2 = taylor2(a,b,n,y0,f,df);
y3 = meuler(a,b,n,y0,f);
y4 = heun(a,b,n,y0,f);
y5 = rk4(a,b,n,y0,f);

fprintf('%5s %11s %11s %11s %11s %11s %11s\n','t_i','Euler','Taylor2','Mod Euler','Heun','RK4','true');
fprintf('-----------------------------------------------------------------------------------\n');
for i=1:n,
 ti=a+i*h;
 fprintf('%5.7g %11.8f %11.8f %11.8f %11.8f %11.8f %11.8f\n',ti,y1(i+1),y2(i+1),y3(i+1),y4(i+1),y5(i+1),y(ti));
end;

% Sample output
%   t_i       Euler     Taylor2   Mod Euler        Heun         RK4        true
% -----------------------------------------------------------------------------------
%   2.5  2.00000000  1.75000000  1.81250000  1.84648277  1.83332336  1.83333333
%     3  2.62500000  2.42578125  2.48155308  2.50941231  2.49997119  2.50000000

