% Find the zero of a function fnc(x) using Bisection
% Animate the search
% [xout,fx,iters] = bisect_anim('fnc',xin,tol)
% 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
%  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)
xlower(1) = xin(1);
xupper(1) = xin(2);
flower(1) = eval([fnc,'(xlower)']); 
fupper(1) = eval([fnc,'(xupper)']);
width = xupper(1)-xlower(1);
xpnt(1)=xlower(1);
fpnt(1)=flower(1);
xpnt(2)=xupper(1);
fpnt(2)=fupper(1);
i=1;
iters = 0;
disp('iteration            interval width')
disp('xlower               xupper')
disp('f(xlower)            f(xupper) = ')
format long e
[[i-1,width];[xlower(i), xupper(i)];[flower(i),fupper(i)]]
pause
i=2;
ageupper(1) = 1;
ageupper(2) = 1;
agelower(1) = 1;
ageupper(2) = 1;
len = 2;
% while (width > tol & len > 1),
while (width > tol),
   xmid = (xlower(i-1) + xupper(i-1))/2;
   fmid = eval([fnc,'(xmid)']);
%  if (fmid*flower(i-1) <= 0),
   if (sign(fmid)*sign(flower(i-1)) <= 0),
      ageupper(i) = i;
      agelower(i) = agelower(i-1);
      xupper(i) = xmid;
      fupper(i) = fmid;
      xlower(i) = xlower(i-1);
      flower(i) = flower(i-1);
   else
      agelower(i) = i;
      ageupper(i) = ageupper(i-1);
      xlower(i) = xmid;
      flower(i) = fmid;
      xupper(i) = xupper(i-1);
      fupper(i) = fupper(i-1);
   end
   width = xupper(i)-xlower(i);
   disp('iteration            interval width')
   disp('xlower               xupper')
   disp('f(xlower)            f(xupper) = ')
   [[i-1,width];[xlower(i), xupper(i)];[flower(i),fupper(i)]]
   pause
   i=i+1;
   iters = iters + 1;
   xpnt(i)=xmid;
   fpnt(i)=fmid;
   if (rem(i,3)==0),
      subs = (max(i-4,1):i);
      subss = (max(i-4,1):i-1);
      plotxmin = min([xpnt(subs),xlower(subss)]);
      plotxmax = max([xpnt(subs),xupper(subss)]);
      w=plotxmax-plotxmin;
      plotxmin = plotxmin - w/20;
      plotxmax = plotxmax + w/20;
      plotxvals = (plotxmin:(plotxmax-plotxmin)/500:plotxmax);
      if (plotxmax==plotxmin), plotxvals = plotxmax; end
      len=length(plotxvals);
      if (len>1), 
         ppts = nnz(plotxvals(2:len)-plotxvals(1:len-1))+1;
      else, 
         ppts = 1;
      end
      plotfvals = eval([fnc,'(plotxvals)']);
      disp(['number of points plotted = ',int2str(ppts)])
%     if (len < 20), [plotxvals',plotfvals'], end
      wf = max([(max(plotfvals)-min(plotfvals))/20, ...
            eps*max(abs(plotfvals)),1e-308]);
      figure(1), hold off, clf
      plot(plotxvals-plotxmin,plotfvals,'g.'), grid
      hold on
      axisarg=[0,plotxmax-plotxmin,min(plotfvals)-wf,max(plotfvals)+wf];
      axisarg(2)=max([axisarg(2),eps*plotxmax,1e-308]);
      axis(axisarg);
      plot(xpnt(subs)-plotxmin,fpnt(subs),'rx')
      plot(xpnt(subs)-plotxmin,fpnt(subs),'ko')
      plot([0,plotxmax-plotxmin],[0,0],'k-')
      for label = subs,
         text(xpnt(label)-plotxmin,fpnt(label)-wf,int2str(label))
      end
      title(['Bisection result for i= ',int2str(min(subs)), ...
             ' to  ',int2str(i)])
      xlabel(['x runs from ',num2str(plotxmin,16),' to ', ...
             num2str(plotxmax,16)])
      pause
   end       
end
xout=[xlower(i-1),xupper(i-1)];
fx=[flower(i-1),fupper(i-1)];