% Find the solutions of the equation f(x) = 0 using Newton iteration method, 
% where f(x) = e^x - x - 2.

clear, clc

% The initial guess
x = 0.5; sol = 1.1461932206205825852;
% x = -0.3; sol = -1.8414056604369606378;

fprintf('initial: x = %18.16f\n',x);
fprintf(['iter.                x                     digits' ...
    '         log(|err|)           order\n']);
dashString = repmat('-', 1, 110);
fprintf('%s\n', dashString);

% Set the tolerance
tol = 1e-15;

% The number of iteration
ite = 0;

% Set the iteration function
g = @(x) x - (exp(x) - x - 2)/(exp(x) - 1);
f = @(x) exp(x) - x - 2;
fprime = @(x) exp(x) - 1;

% The plotting part
LW = 'LineWidth'; lw1 = 2; lw2 = 0.5;
LS = 'LineStyle'; ls = '--'; 
FS = 'fontsize'; fs = 12;
figure;
set(gcf,'position',[600,85,900,650])
fplot(f, LW ,lw1); hold on  % Plot the function
xmin = min(x - 0.5, -3); xmax = max(x + 0.5, 2);  % Get the region
axis([xmin xmax -2.2 f(xmax)]);
xscale = xmax - xmin; yscale = f(xmax) + 2.2;
title('Newton')
scatter(x,0,'filled','r'); hold on
xite = ['x' sprintf('%d',ite)];
text(x-0.01*xscale, -yscale*0.022, xite, FS, fs);
plot([xmin, xmax],[0, 0], LW, lw2, 'Color', 'black'); hold on

err = 0; % store the sequence of errors
% do the iteration
while 1
    pause
    ite = ite + 1;
    xnew = g(x);
    err = [err;abs(xnew - x)];
    logerr = log10(err(ite + 1));
    digits = max(0, - floor(log10(abs(xnew - x))) - 1);
    if ite > 2
        order = log(err(ite + 1)/err(ite))/log(err(ite)/err(ite - 1));
        if order == Inf
            fprintf('%2d   %20.16f  %5.2d  %18.8f    \n',...
                ite, xnew, digits, logerr);
        else
            fprintf('%2d   %20.16f  %5.2d  %18.8f     %8.3f\n',...
                ite, xnew, digits, logerr, order);
        end
    else
        fprintf('%2d   %20.16f    %25.8f     \n',...
            ite, xnew,logerr);
    end
    plot([x, x],[0, f(x)], LW, lw2, LS, ls, 'Color', 'black'); hold on
    fplot(@(t) fprime(x)*(t - x) + f(x), LW, lw2, 'Color', 'black'); hold on
    scatter(xnew,0,'filled','r'); hold on
    if abs(xnew - x) > 0.03*xscale
        xite = ['x' sprintf('%d',ite)];
        text(xnew-0.01*xscale, -yscale*0.022, xite, FS, fs);
    end
    if abs(xnew - sol) < tol
        fprintf('%s\n', dashString);
        fprintf('True solution: x = %18.16f\n',sol)
        break
    else
        x = xnew;
    end
end