% Root finding using BISECTION METHOD. 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.

 clear;  
 format short e;
 fprintf(' Bisection Method.\n');
 
	disp(' Introduce the ends A and B. Choose them such that A < B');
   A=input('A=');
   B=input('B=');     
   fA=f(A);
   SA=sign(fA);
   fB=f(B);
   SB=sign(fB);
   
   while SA==SB
   	disp('There is not root in the given interval');
   	disp(' Choose new values for the ends of the interval');
   	A=input('A=');
   	B=input('B=');
    fA=f(A);
    SA=sign(fA);
    fB=f(B);
    SB=sign(fB);	
   end; 

% Ends of the original interval where the root is located
   A0=A;
   B0=B;
   % Input a tolerance for the error "TOL" and a maximum number of
   % iterations "itmax".
   TOL=input('TOL=');
   itmax=input('itmax=');

   for it=1:itmax	
		  p=(A+B)/2;
          fp=f(p);
          Sp=sign(fp);
          Coterr=(B-A)/2;
          disp('        it         A             B         Cterr          p            f(p)');
				disp( [it, A, B, Coterr, p, f(p)]);
   			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         Cterr          p            f(p)');
				disp( [it, A, B, Coterr, p, fp]);

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

x1=-5:.05:3;
y1=f(x1);
plot(x1,y1)
xlabel('x'); ylabel('f(x)');
grid on

% Exact solution

ExValue=fzero('f', 1);
disp('Exact Value');
format long;
disp(ExValue);