%%Bisection Method for Root Finding

% Idea: Start with two points, xleft and xright with different signs
% (so there is a root between them!).
% Take the midpoint of the interval and choose which half of the
% interval to keep.

% Specifically: Start with xleft and xright where f(xleft) and
% f(xright) have different signs. Then, at each iteration
% Set xnew = (xleft + xright)/2;
% if sign(f(xleft)) == sign(f(xnew))
%   xleft = xnew;
% else
%   xright = xnew;
% end

%Repeat this process until the estimated error is below epsilon
%(our desired accuracy of the solution).
%We halve the error at each iteration, since our interval shrinks
%by a factor of 1/2.
%This means that we have a linear rate of convergence!

%Example: f(x) = 6 - 2x^2
%True roots are at +-sqrt(3)
f = @(x)6 - 2*x.^2;

%First step: plot to find a good interval
xx = -3 : 0.01 : 3;
yy = f(xx);
ymin = -7;
ymax = 7;
xmin = -3;
xmax = 3;
plot(xx,yy,'b',[0,0],[ymin,ymax],'k',[xmin,xmax],[0,0],'k')
%Adds a grid to help see exactly where the line is
grid on
axis([xmin xmax ymin ymax])

% Set up to call bisection
xleft = 1;
xright = 2;
epsilon = 1e-12;  % desired accuracy
root = bisection(f,xleft,xright,epsilon);

% How accurate were we?
precision = 20; % precision to use in num2str in this script
disp(['Actual root = ', num2str(sqrt(3), precision)])
disp(['Calculated root = ', num2str(root, precision)])

