function bisect

%  Solve  f(x) = 0  using the bisection method
%  Usese exact solution to check on error

%  Input:
%	a = left endpoint of interval
%	b = right endpoint of interval
%	tol = tolerance for stopping bisection method
%	f(x) This is at end of file

a=1/5;
b=3;
tol=10^(-6);

% exact solution
xe=2;

% plot function for a < x < b to show there is a solution
figure(1)
clf
x=linspace(a,b,200);
for ix=1:200
    y(ix)=f(x(ix));
end
clf
set(gcf,'Position', [1 925 560 420])
plot(x,y,'r','LineWidth',1.2)
hold on
plot([a b],[0 0],'k','LineWidth',1.2)
xlabel('x-axis')
ylabel('y-axis')
grid on
box on
set(gca,'FontSize',16,'FontWeight','bold')
pause

% bisection method
mid=0.5*(a+b);
it=0;
errM=0.5*(b-a);
while (b-a)>tol
    fmid=f(mid);
    if f(a)*fmid<0
        b=mid;
    elseif f(b)*fmid<0
        a=mid;
    else
        break
    end 
    mid=0.5*(a+b);    
    it=it+1; iteration(it)=it; error(it)=abs(xe-mid);
    fprintf('\n %d  Computed Solution = %14.8e    Error = %5.2e',it,mid,error(it))
    pause
end
fprintf('\n\n')

% plot error curve
figure(2)
clf
set(gcf,'Position', [3 431 560 420])
semilogy(iteration,error,'--or','LineWidth',1.5,'MarkerSize',8)
xlabel('Iteration Step (n)')
ylabel('Error')
grid on
box on
set(gca,'FontSize',16,'FontWeight','bold')



function g=f(x)
g=x*(x-2)*(x-4);