function [r,iter] = bisection(f,a,b,tol,nmax)
  % by Shelvean Kapita
  % bisection algorithm for determining roots
  % example usage
  % f = @(x) x - cos(x);
  % a = 0;
  % b = 1;
  % tol = 1e-8;
  % nmax = 100;
  % [r,iter] = bisection(f,a,b,tol,nmax)

  iter = 1; % initialize count
  r = zeros(nmax,1); % initialize vector of iterates
  while ((b-a)/2>tol) && (iter < nmax)
      c = (b+a)/2; % midpoint of [a,b]
      r(iter) = c; % update iterate
      if f(c)==0 % c is the solution!
          break
      elseif f(a)*f(c)<0
          b = c;
      else
          a = c;
      end
      iter = iter + 1;
  end
  r = r(1:iter-1); % return iterates
  iter = iter-1;