% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 16: Bisection method for solving f(x) = 0 in [a,b]
% A. Meseguer (2007)

% On input
% 1. [a,b]: interval (it asumes that f(a)f(b) < 0)
% 2. tol:   tolerance so that abs(bk-ak)/2 < tol
% 3. itmax: stop criteria for exceeding number of iterations
% 4. fun:   function to be solved


% On output
% 1. xk:   resulting iteration points
% 2. resd: resulting residuals iteration
% 3. it:   number of required iterations to satisfy tolerance 

  function [xk,resd,it] = bisection(a,b,tol,itmax,rootfun)

  ak = a ; bk = b ; xmid = (ak + bk)*.5 ; it = 0 ;
  xk = [xmid] ; fmid = feval(rootfun ,xmid) ; resd = [abs(fmid)] ;
  tolk = .5*abs(bk-ak) ;

  while it < itmax & tolk > tol
   fa = feval(rootfun,ak) ;

   if fmid*fa < 0
    bk = xmid ;
   else
    ak = xmid ;
   end
   tolk = .5*abs(bk-ak) ;
   xmid = (ak + bk)*.5 ; fmid = feval(rootfun ,xmid) ;
   xk = [xk xmid] ; resd = [resd abs(fmid)] ; it = it + 1 ;
  end





