% bisectiondemo
% Demonstration of Bisection Root Finding Algorithm
%
% Usage: Write the function to be solved in a separate .m file
% e.g. f42.m
% y = x.^3 + x.^2 - 10;
% 
% Then call bisectiondemo2 to find its root within the interval [1,2] 
%  and upper bound for E
clear all
close all
clc

% For this example, we want to find the actual root between [1,2] to
%  calculate the true error

r = roots([1, 1, 0, -10]);  %Matlab's function to find roots of polynomials
Vt = r(3);                  % True value = 1.8675...

a = 1;
b = 2;

% First plot to see the behaviour
fplot(@(x) f42(x),[a,b],'LineWidth',2)
hold on
xlabel('x')
ylabel('f(x)')
title('f(x) = x^3 + x^2 - 10')
fplot(0,[a,b],'-k')
grid on


% initialize other variables
maxit = 500;    % maximum number of iterations
fa = f42(a);
fb = f42(b);
es = 0.0005;

% test to ensure a single root is bracketed
if (sign(fa)==sign(fb))
    fprintf('\n Error: no root or even number of roots in interval.\n');
    return;     % terminates the function
end

% prepare the output header
fprintf('\n  #         xi          a          b       er       Actual E       Upper Bound for E\n\n');

count = 0;  % iteration counter
er = 1;     % initial error value to enter loop
xi = a;     % initialize xi and xold
xold = xi;

% % main loop
while (abs(er) > es && count<maxit)
    count = count +1;
    xold = xi;  % keep the previous estimate
    
    xi = (a + b)/2; % calculate the new estimate
    fxi = f42(xi);
    plot(xi,fxi,'o','MarkerEdge',[0.6350, 0.0780, 0.1840]);  % allows us to visualize the approximations
       
    if count ~= 1
        er = (xi-xold)/xi;  % relative error
    end
    Et = Vt - xi;
    L = b - a;
    
    fprintf('%3g %10g %10g %10g %10.4g %12.4g %15.4g\n',count-1,xi,a,b,er,Et,L);
        
    product = fa*fxi;   % test product
    
    if product == 0     % root has been found
        er = 0;
    elseif product < 0  % root in fist half
        b = xi;
        fb = fxi;
    else                % root in second half
        a = xi;
        fa = fxi;
    end
end

disp(' ')
disp(['Root has been found at x = ',num2str(xi)])
disp(' ')
disp('--------------------------------------------')
disp(' ')


