% Initial value
y0 = 0.1;

% Initial time and end final time
t0 = 0;
T = 10;

% Number of time steps
nt = 30;

% Compute discretization width
h = (T - t0) / nt;

f = @(t,y) y.*(1-y);

t = linspace(t0,T,nt+1);

% Compute true solution
y_true = 1 ./ (1 + (1/y0-1) * exp(-t));

% Compute numerical solution
y_euler = euler_explicit(h,f,t0,y0,T);
% y_mod = euler_modified(h,f,t0,y0,T);
y_improved = euler_improved(h,f,t0,y0,T);

figure(1), clf, hold on
plot(t,y_true,'-k','DisplayName','True Solution')
plot(t,y_euler,'-b','DisplayName','Explicit Euler')
% plot(t,y_mod,'-r' ,'DisplayName','Modified Euler')
plot(t,y_improved,'-r' ,'DisplayName','Improved Euler')
legend('Location','SouthEast')

err = zeros(2,10);

for i = 1:10,
	h = 2^(-i);
	t = t0:h:T;
	y_true = 1 ./ (1 + (1/y0-1) * exp(-t));
	y_euler = euler_explicit(h,f,t0,y0,T);
	% y_mod = euler_modified(h,f,t0,y0,T);
	y_improved = euler_improved(h,f,t0,y0,T);
	H(i) = h;
	err(1,i) = norm( y_euler - y_true, Inf);
	% err(2,i) = norm( y_mod - y_true, Inf);
	err(2,i) = norm( y_improved - y_true, Inf);
end

figure(2), clf
loglog(H, err(1,:), '-b','DisplayName','Explicit Euler')
hold on
loglog(H, err(2,:), '-r','DisplayName','Modified Euler')
axis equal
legend('Location','SouthEast')

fprintf('Estimated order of convergence  for explicit Euler: %5.4f\n', eoc(H', err(1,:)'));
fprintf('Estimated order of convergence  for modified Euler: %5.4f\n', eoc(H', err(2,:)'));
