% Root finding using BISECTION METHOD as a SUBROUTINE. 
% This code is an
% improvement of the Bisection code. We have incorporated 
% a change to evaluate the functions only once at the needed
% points. Also the function sign(x) has been used to avoid 
% possible overflow or underflow when working with too small 
% or too large numbers.

function Bisection(A,B,TOL,itmax)
 format short e;
 fprintf(' Bisection Method.\n'); 
   fA=f(A);
   SA=sign(fA);
   fB=f(B);
   SB=sign(fB);
   
   if SA==SB
   	disp('There is not root in the given interval');
   	disp(' Choose new values for the ends of the interval');
    break;
   end; 

% Exact solution

ExValue=fzero('f', 10^(-6));
disp('Exact Value');
disp(ExValue);
   
% Ends of the original interval where the root is located
   A0=A;
   B0=B;
   for it=1:itmax	
		  p=(A+B)/2;
          fp=f(p);
          Sp=sign(fp);
          Coterr=(B-A)/2;
          Trueerror=abs(ExValue-p);
          disp('        it         A             B         errbound          p          f(p)       AbsErr');
				disp( [it, A, B, Coterr, p, f(p),Trueerror]);
   			if fA==0, fprintf('approx. root= %12.5f\n',A), break, end;
   			if fB==0, fprintf ('approx. root=%12.5f\n', B), break, end
            if fp==0 | Coterr <TOL,      		
               fprintf('approx. root= %12.5f\n',p), break, end
   			if SA==Sp, 
                A=p;
                fA=fp;
                SA=Sp;
   			else 
                B=p;
                fB=fp;
                SB=Sp;
   			end       
end
disp('.................................................................................');        
disp('Final Results');
disp('        it         A             B         errbound          p          f(p)');
				disp( [it, A, B, Coterr, p, fp]);

% Function defining the equation for which roots are being found.

x1=-4:0.01:3;
y1=f(x1);
plot(x1,y1)
xlabel('x'); ylabel('f(x)');
grid on