% sample function to find the zero of a function using Bisection
% [xout,fx,iters] = bisect('fnc',xin,tol,debug)
% Inputs:
%    fnc - name of function
%    xin(1:2) - lower and upper limits of interval to search for a zero
%                func must have different signs at xin(1),xin(2)
%    tol - stop refining interval when its width is < tol
%    debug - print results of each step when debug>0
%  Outputs:
%    xout(1:2) - output interval "guaranteed" to contain a zero
%    fx(1:2) - value of func at xout(1),xout(2)
%    iters = number of iterations to converge
function [xout,fx,iters] = bisect(fnc,xin,tol,debug)
xlower= xin(1);
xupper = xin(2);
flower = eval([fnc,'(xlower)']); 
fupper = eval([fnc,'(xupper)']);
width = xupper-xlower;
iters=0;
if (debug>0),
   disp('iteration            interval width')
   disp('xlower               xupper')
   disp('f(xlower)            f(xupper) = ')
   format long e
   [[iters,width];[xlower,xupper];[flower,fupper]]
   disp('pause - hit return to continue'), pause
end
while (width > tol),
   xmid = (xlower + xupper)/2;
   fmid = eval([fnc,'(xmid)']);
   if (fmid*flower <= 0),
      xupper = xmid;
      fupper = fmid;
   else
      xlower = xmid;
      flower = fmid;
   end
   width = xupper-xlower;
   iters=iters+1;
   if (debug>0),
      disp('iteration            interval width')
      disp('xlower               xupper')
      disp('f(xlower)            f(xupper) = ')
      format long e
      [[iters,width];[xlower,xupper];[flower,fupper]]
      disp('pause - hit return to continue'), pause
   end
end
xout=[xlower,xupper];
fx=[flower,fupper];