%  Demo_bisection.m
%
%  Solve F(x) = exp(x) - 5 = 0 using bisection method.

  disp(' Use bisection method to solve F(x) = exp(x) - 5 = 0.')
  format long;      %  Display 15 digits instead of usual 5.
  x_true = log(5)
  %  initial interval [a0,b0] should contain true solution.
  a0 = 0;
  b0 = 3;
  xmid = (a0 + b0) / 2;
  maxiter = 15; %%%input(' Max. no. of Newton iterations = ');
  Ftol = 1e-6; %%%input(' Absolute stopping tolerance, Ftol = ');

%  Save solutions and errors in xvec and evec.

  xvec = [];
  evec = [];

%  Plot function and initial left and right endpoints.
  
  t = [a0:(b0-a0)/100:b0]';
  F = exp(t) - 5;
  figure(1)
    clf     %  Clear previous graph
    plot(t,F,'-', a0,0,'o', b0,0,'*')
    xlabel('x axis')
    legend('F(x)', 'left endpoint', 'right endpoint')
  hold on  %  Hold graph; subsequent plot commands plot over it instead
           %  of resetting it.

%  Initialization.

  k = 1;                                   % Iteration count.
  ak = a0;                                 % Left endpoint
  bk = b0;                                 % Right endpoint
  Fak = exp(ak) - 5;                       % F(ak)
  Fbk = exp(bk) - 5;                       % F(bk)
   
%  iteration

  while k < maxiter & abs(Fak) > Ftol,
    
    %  Compute new approximation x_new to solution and plot line segment 
    %  between (xk,Fk) and (x_new,0).
    
    xmid = (ak+bk) / 2;
    Fxmid = exp(xmid) - 5;
    if Fxmid*Fak < 0   % F(ak) and F(xmid) have opposite signs
      bk = xmid;       % Reset right endpoint
      Fbk = Fxmid;
    else               % F(ak) and F(xmid) have same signs
      ak = xmid;       % Reset left endpoint
      Fak = Fxmid;
    end

    xk = xmid;
    k = k + 1;
    
    %  Display results.

    plot(t,F,'--', ak,0,'o', bk,0,'*')
    xlabel('x axis')
    legend('F(x)', 'left endpoint', 'right endpoint')
    fprintf(' Hit any key to continue.\n');
    pause

    %  Save approximate solution and error.
    
    xvec(k) = xk;
    evec(k) = abs(xk - x_true);
  end

%  Display vectors of approximate solutions and errors.

  fprintf('\n Iter   Approx       Error\n')
  k_last = k;
  for k = 1:k_last
    fprintf(' %2.0f   %10.8f   %6.4e\n', k, xvec(k), evec(k));
  end

  figure(2)
    semilogy(evec,'o-')
    xlabel('Bisection Iteration')
    title('Absolute error of approximation vs bisection iteration')
    grid