% Root finding using BISECTION METHOD.

 clear;  
 format short e;
 fprintf(' Bisection Method.\n');

 % *****************************DATA INPUT*********************************
 
    disp(' Introduce the ends A and B. Choose them such that A < B');
    A=input('A=');
    B=input('B=');
    Product=f(A)*f(B);
    while Product>0
   	    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=');
   	    Product=f(A)*f(B);
	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=');
   
% *************************************MAIN LOOP***************************

   for it=1:itmax	
		  p=(A+B)/2;
          Coterr=(B-A)/2;
          disp('        it         A             B         Cterr          p            f(p)');
				disp( [it, A, B, Coterr, p, f(p)]);
   			if f(A)==0 
                fprintf('exact root= %f   Number of iterations=%d\n',A,it), break
            end;
   			if f(B)==0
                fprintf('exact root= %f   Number of iterations=%d\n',B,it), break
            end
            if f(p)==0 | Coterr <TOL,
                fprintf('approx. root= %f   Number of Iterations=%d\n',p,it), break
            end
   			if f(A)*f(p)>0
                A=p;
   			else 
                B=p;
   			end       
   end

% *************************************FINAL OUTPUT************************

disp('.................................................................................');        
disp('Final Results');
disp('        it         A             B         Cterr          p            f(p)');
disp( [it, A, B, Coterr, p, f(p)]);

% Graph of a function that defines the equation for which roots are being found.

step=(B0-A0)/100;
x1=A0:step:B0;
y1=f(x1);
plot(x1,y1)
xlabel('x'); ylabel('f(x)');
grid on

% Exact solution. MATLAB root finding command "fzero".

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