function [left_point right_point] = bisection(f,a,b,err)
% [left_point right_point] = bisection(f,a,b,err)
% applies the bisection method for the solution of the function f(x) = 0 on
% the interval [a,b] with accuracy err.
%
% input:
% f   ...function handle to f
% a   ...left endpoint of the initial interval
% b   ...right endpoint of the initial interval
% err ...desired accuracy
%
% The method requires that f(a) and f(b) have a different sign.
%
% Example:
%
% % compute the solution of tan(x) = 0 on [2,4] to an accuracy of 10^-12
% [a b] = bisection(@tan,2,4,10^(-12))
% % as expected the distance between the solution and pi is smaller than
% the error 10^-12:
% abs(b-pi), abs(a-pi)
%
% % Compute a solution of 3cos(x)sin(x^2) = 1 on [0,10].
% % Note that this equation has several solutions on this interval and it
% % is not clear which solution the method selects.
% fun = @(x)(3*cos(x)*sin(x^2)-1);
% [a b] = bisection(fun,0,10,10^(-12))


% check whether [a,b] might be an interval
if ~(a < b)
    error('The input points have to form a non-trivial interval.')
end

% compute f at the boundaries of the interval
f_left = f(a);
f_right = f(b);

% compute the current error
err_current = b-a;

% check whether the input already solves the equation
if f_left == 0
    left_point = a;
    right_point = a;
    flag_converged = 1;
elseif f_right == 0
    left_point = b;
    right_point = b;
    flag_converged = 1;    
else
    % initialise the interval
    left_point = a;
    right_point = b;
    sgn_left = sign(f_left);
    sgn_right = sign(f_right);
    % check whether the signs of f on the boundaries are different
    if sgn_left == sgn_right
        error('The signs of the function values at the endpoints of the interval need to be different.');
    end
    flag_converged = 0;
end

% main iteration
while ~flag_converged
    % compute the midpoint of the current interval
    midpoint = (left_point+right_point)/2;
    % evaluate f at the midpoint
    f_mid = f(midpoint);
    % check whether the midpoint is a solution
    if f_mid == 0
        left_point = midpoint;
        right_point = midpoint;
        flag_converged = 1;
    % if the sign of the function value at the midpoint is the same as at
    % the left endpoint, then replace the left endpoint by the midpoint.
    elseif sign(f_mid) == sgn_left
        left_point = midpoint;
        % update the length of the interval
        err_current = err_current/2;
        % stop the iteration if the desired accuracy is reached
        flag_converged = (err_current < err);
    % at that stage, the function value at the midpoint should have the
    % same sign as the value at the right endpoint. Thus the right endpoint
    % should be replaced by the midpoint.
    else
        right_point = midpoint;
        % update the length of the interval
        err_current = err_current/2;
        % stop the iteration if the desired accuracy is reached
        flag_converged = (err_current < err);
    end        
end