function [r,sequence_x,sequence_fx]=bisection(f,a,b,n_max,bisection_tolerance)
% function [r,sequence_x,sequence_fx]=bisection(f,a,b,n_max,bisection_tolerance)
% runs the bisection method starting with two points, 
% a and b, to find a zero r of the function f with the tolerance 
% bisection_tolerance in n_max steps, whichever is reached first. 
% It returns the zero r, and the sequence of approximations 
% sequence_x as well as the function values at these points sequence_fx
% The first two values of sequence_x are a and b. 

%Revision 1.1 by Andrew Knyazev. Oct. 14, 2002. 

sequence_x=zeros(n_max+2,1); sequence_fx=zeros(n_max+2,1); %to avoid dynamic memory allocation

fa=feval(f,a);fb=feval(f,b);
fprintf('Initial values: a= %e , b= %e , f(a)= %e ,f(b)= %e \n \n',a,b,fa,fb);

sequence_x(1)=a; sequence_x(2)=b;
sequence_fx(1)=fa; sequence_fx(2)=fb;

if sign(fa)==sign(fb) 
    error('The function has the same sign at a and b. Abort bisection.');
end

bisection_error=b-a;

for n=1:n_max
    bisection_error=bisection_error/2;
    c=a+bisection_error;
    fc=feval(f,c);
    sequence_x(n+2)=c; sequence_fx(n+2)=fc;
    fprintf('Iteration n= %i , c= %e , f(c)= %e ,bisection_error= %e \n',n,c,fc,bisection_error);
    if abs(bisection_error) < bisection_tolerance
        fprintf('Converged within the tolerance %e . The zero found is %e \n',bisection_tolerance,c);
        r=c;
        sequence_x(n+2:n_max+2)=[]; sequence_fx(n+2:n_max+2)=[]; %remove artificial zeros
    return
    end
    if sign(fa)~=sign(fc)
        b=c; fb=fc;
    else
        a=c; fa=fc;        
    end
end
r=c;
return