% bisectiondemo
% Demonstration of Bisection Root Finding Algorithm
%
% Usage: Write the function to be solved in a separate .m file
% e.g. f41.m or f42.m
% 
% Then call bisectiondemo to find its roots
clear all
close all
clc

%% First Root
a = -0.6;
b = -0.5;

% First plot to see the behaviour
fplot(@(x) f41(x),[-1,1],'LineWidth',2)
hold on
xlabel('x')
ylabel('f(x)')
title('f(x) = sin(5x)+ cos(2x)')
fplot(0,[-1,1],'-k')
grid on


% initialize other variables
maxit = 500;    % maximum number of iterations
fa = f41(a);
fb = f41(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  #         a          b          xi       f(a)       f(b)       f(xi)        er\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 = f41(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
    
    fprintf('%3g %10g %10g %10g %10.3g %10.3g %10.3g %10.4g\n',count-1,a,b,xi,fa,fb,fxi,er);
        
    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(' ')


%% Second Root
a = -0.3;
b = -0.2;

% initialize other variables
fa = f41(a);
fb = f41(b);

% 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  #         a          b          xi       f(a)       f(b)       f(xi)        er\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 = f41(xi);
    plot(xi,fxi,'o','MarkerEdge',[0.8500, 0.3250, 0.0980]);  % allows us to visualize the approximations
  
    
    if count ~= 1
        er = (xi-xold)/xi;  % relative error
    end
    
    fprintf('%3g %10g %10g %10g %10.3g %10.3g %10.3g %10.4g\n',count-1,a,b,xi,fa,fb,fxi,er);
    
    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(' ')


%% Third Root
a = 0.6;
b = 0.7;

% initialize other variables
fa = f41(a);
fb = f41(b);

% 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  #         a          b          xi       f(a)       f(b)       f(xi)        er\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 = f41(xi);
    plot(xi,fxi,'o','MarkerEdge',[0.4940, 0.1840, 0.5560]);  % allows us to visualize the approximations
    
    
    if count ~= 1
        er = (xi-xold)/xi;  % relative error
    end
    
    fprintf('%3g %10g %10g %10g %10.3g %10.3g %10.3g %10.4g\n',count-1,a,b,xi,fa,fb,fxi,er);
    
    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(' ')
