function p = bisect(f,a,b,tol,nmax);
%        p = bisect(f,a,b,tol,nmax);
%  Performs the bisection method to approximate a zero of a function f(x).
%
%  Input: f (string) name of the function 
%         a, b   endpoints of the initial interval
%         tol    iteration stops when (b-a)/2 < tol
%         nmax   maximum number of permitted iterations
%  Output: p the final approximation of a root (the midpoint of the
%              final interval)

%  A table of results for each iteration is also displayed.
%

   if (a >= b), error(' a must be less than b'); end;

%  Note the use of the MATLAB command `feval' to evaluate a function 
%  whose name is passed in as an input parameter.
   fa = feval(f,a); fb = feval(f,b);

   if (fa*fb >= 0), 
      error(' invalid [a,b]; f(a) and f(b) do not have opposite signs');
   end;

   n=1;
   disp(sprintf('%4s  %12s  %12s %12s %12s %12s','n', 'a ','b', 'p', 'f(p)', '(b-a)/2'));
   while (n <= nmax),
      p = a + (b-a)/2;
      fp = feval(f,p);
      
      disp(sprintf('%4d  %e   %e  %e  %e  %e',n,a,b,p,fp,(b-a)/2));

      if (fp == 0 | (b-a)/2 < tol), return; end;
      n=n+1;
      if (fa*fp > 0),
          a = p; fa = fp;
      else
          b=p;
      end;
   end;
   if (n > nmax), disp('did not congerge in nmax iterations'); end;