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

clear, clc

% The initial guess
% x0 = -1; x1 = -0.2; sol = -1.8414056604369606378; 
x0 = 0.5; x1 = 3; sol = 1.1461932206205825852; 

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

% Set the tolerance
tol = 1e-14;

% The number of iteration
ite = 0;

% Set the function
f = @(x) exp(x) - x - 2;
fold = f(x0);

% 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([x0 - 0.5, x1 - 0.5, -3]); xmax = max([x0 + 0.5, x1 + 0.5, 2]); 
axis([xmin xmax -2.2 f(xmax)]);
xscale = xmax - xmin; yscale = f(xmax) + 2.2;
title('Secant')
plot([xmin, xmax],[0, 0], LW, lw2, 'Color', 'black'); hold on
scatter([x0, x1],[0, 0],'filled','r'); hold on
text(x0-0.01*xscale, -yscale*0.022, 'x0', FS, fs);
text(x1-0.01*xscale, -yscale*0.022, 'x1', FS, fs);

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