% Find a root of any function
% Function specified in fun.m and evaluated by fun(x) where x is any real value

lower_limit = input('What is the lower limit for searching? '); % Get user input for lower limit
upper_limit = input('What is the upper limit for searching? '); % Get user input for upper limit
lower_value = fun (lower_limit); % Evaluate the function at the lower limit
upper_value = fun (upper_limit); % Evaluate the function at the upper limit
if upper_value == 0 % Check to see if upper limit lies on root
    display ('Found Root at ') % Then found root
    rootat = upper_limit % Display value of root
    break
end
if lower_value == 0 % Check to see if lower_limit lies on root
    display ('Found root at ') % Then found root at
    rootat = lower_limit % Display value of root
    break
end

if (lower_value < 0 & upper_value < 0) | (lower_value > 0 & upper_value > 0) % Check to see if a root is obvious between given limits
    display ('There does not appear to be a root') % Will change this so that a different algorithm will be used
    break
end

for i = [1 : 1 : 1000] % Main Loop.  Set up 1000 iterations, should be enough.  If no root found then will resort to another algorithm
    middle_limit = (lower_limit + upper_limit)/2; % Evaluate the middle limit - pure bisection
    middle_value = fun (middle_limit); % Evaluate the value at this middle limit
    
    % Lets check to see if we have landed on a root and display that one if we have
    if middle_value == 0 % Check to see if middle_limit lies on root
        display ('Found root at ') % Then found root at
        rootat = middle_limit % Display value of root
        break
    end    

    % Lets find where the root lies and adjust limits accordingly
    if (lower_value < 0 & middle_value > 0) | (lower_value > 0 & middle_value < 0) % If root lies between lower and middle limits
        upper_limit = middle_limit; % then discard upper limit
        upper_value = fun (upper_limit); % Evaluate the function at the new upper limit
    end
    if (middle_value < 0 & upper_value > 0) | (middle_value > 0 & upper_value < 0) % If root lies between upper and middle limits
        lower_limit = middle_limit; % then discard lower limit
        lower_value = fun (lower_limit); % Evaluate the function at the new lower limit
    end
    
    % Lets check to see if we are close enough to the root
    closeness = abs(lower_value) + abs(upper_value); % Check to see how close the values are to zero
    if closeness < 0.001 % IF they are within threshold
        display ('Found Root at ') % Then found root
        rootat = (lower_limit + upper_limit) / 2 % Display value of root
        break
    end
end