% WARNING!  This is for f(x) = (x+1)(x-1)^2.  If you're going to
% change the function, you have to change it in _three_ places!
%
% [root,value,steps,A] = bisect(a_0,b_0,tol) reads in a left-hand endpoint
% a_0, a right-hand endpoint b_0, and a desired tolerance tol.  it
% bisects until it gets |f(root)| < tol.  It then returns the root,
% the value of the function at the root, and the number of steps it
% took to find it

% A is a diagnostic matrix.  At the moment, it's 100x3 and contains
% A(i,1) = a_i, A(i,2) = c_i, A(i,3) = b_i

function [root,value,steps,A] = bisect(a_0,b_0,tol)

f_test = 2*tol;
steps = 0;

A = zeros(100,3);

while ( abs(f_test) > tol) & (steps < 100) 

% print the interval length to the screen:  (the semi-colon's missing, that's
% why it prints to the screen)
    b_0 - a_0

% change the function here:
    f_left = (a_0+1)*(a_0-1)^2;
% change the function here:
    f_right = (b_0+1)*(b_0-1)^2;

    c_0 = (a_0 + b_0)/2;

% update the counter and fill the diagnostic matrix
    steps = steps + 1;
    A(steps,1) = a_0;
    A(steps,2) = c_0;
    A(steps,3) = b_0;

% change the function here:
    f_test = (c_0 + 1)*(c_0-1)^2;

    if f_left*f_test <= 0
% then there must be a zero between a_0 and c_0, so take c_0 to be 
% the new right-hand endpoint
       b_0 = c_0;
    else
% then there must be a zero between c_0 and b_0, so take c_0 to be 
% the new left-hand endpoint     
       a_0 = c_0;
    end


end


root = c_0;
value = f_test;