% secant
function xr = secant(fun,x0,x1)
% Demonstration of Secant Method
%
% The input function should be of the form f(x)
%       
% Then call secant to find its roots with two initial guesses, x0 and x1
% xr = secant(@fun,x0,x1)
%
% For Example 4.9 use:
% secant(@(x) sin(5.*x)+cos(2.*x), @(x), 0.5, 0.4) 
%


close all
clc

pmin = 0; % x limits of the plot. Change for different functions if needed.
pmax = 1;

% First plot f(x)to see the behaviour
figure
hold on
fplot(@(x) fun(x),[pmin,pmax],'LineWidth',2)
fplot(0,[pmin,pmax],'-k')
ylim([-1.5,2]) % y limits of the plot. Change for different functions if needed.
xlabel('x')
ylabel('f(x)')
title('Secant Demo')
grid on

% Initialize variables
es = 0.0005;
maxit = 100;

% prepare the output header
fprintf('\n  #         xi          er\n');

count = 0;  % iteration counter
xi_1 = x0;  % initialize the estimates
xi = x1;

er = 1;  % relative error

fprintf('%3g %10.6g %10.4g\n',count,xi_1,er);
fprintf('%3g %10.6g %10.4g\n',count+1,xi,er);

h1 = plot([xi_1,xi],[0,0],'o','MarkerEdge',[0.8500, 0.3250, 0.0980]);  % allows us to visualize the approximations
erplot = zeros(1,5);
pause % wait for keypress

count = count+1; % update counter so the first two guesses are iterations 0 and 1

% main loop
while (abs(er) > es && count<maxit)
    count = count +1;
    xi_2 = xi_1; % discard the oldest estimate and update x_(i-1)
    xi_1 = xi;   

    fxi_2 = fun(xi_2); % evaluate the function at the estimates
    fxi_1 = fun(xi_1);
        
    xi = xi_1 - fxi_1*(xi_1 - xi_2)/(fxi_1 - fxi_2); % calculate the new estimate
    
    delete(h1); % for plotting. Removes the old lines on the plot
    if(count>2)
        delete(h2);
        delete(h3);
    end
    h1 = plot([xi_2,xi_2,xi_1],[0,fxi_2,fxi_1],'--','Color',[0.8500, 0.3250, 0.0980]);  % allows us to visualize the approximations
    h2 = plot([xi_1,xi_1],[0,fxi_1],'--','Color',[0.8500, 0.3250, 0.0980]);
    h3 = plot([xi_1,xi],[fxi_1,0],'--','Color',[0.8500, 0.3250, 0.0980]);
       
    er = (xi - xi_1)/xi;  % relative error
    erplot(count) = er;
   
    fprintf('%3g %10.6g %10.4g\n',count,xi,er);
    pause  % Wait for key press
end

if count >= maxit
    disp(' ')
    disp(['Failed to converge in ',num2str(maxit),' iterations'])
    disp(' ')
    disp('--------------------------------------------')
    disp(' ')
    xr = NaN;
else
    disp(' ')
    disp(['Root has been found at x = ',num2str(xi)])
    disp(' ')
    disp('--------------------------------------------')
    disp(' ')
    xr = xi;
end

figure
hold on
fplot(0,[0,count],'-k')
plot([2:count],erplot(2:count),'LineWidth',2)
title('Relative Error')
xlabel('Iteration Number')
ylabel('Relative Error (er)')
hold off

