function x = bisection_root (f, x0, x1, iter)

  % f:        a given scalar function y=f(x)
  % x0:       first starting point
  % x1:       second starting point
  % iter:     stop after maxiter iterations

  left   = x0;
  right  = x1;
  fleft  = feval (f,left);
  fright = feval (f,right);
  if fleft*fright > 0
    error ('Interval does not bracket a root.');
  end

  for j=1:iter
    mid  = (left+right)/2;
    fmid = feval (f,mid);
    if fleft*fmid > 0
      left = mid; 
    else
      right = mid; 
    end
  end
  x = (left+right)/2;

end

  