function [ xnew ] = bisection(func,xleft,xright,epsilon)
%BISECTION Finds roots of func using the bisection method
% xleft and xright define the right and left side of the initial interval
% finds a root within the desired accuracy epsilon
  done = 0;
  iteration = 0;
  
  fprintf('iter %10s  %10s  %10s  %10s  %10s  %10s\n', ... 
'x_left','f(x_left)','x_right','f(x_right)','x_new','f(x_new)')
  
  %Throws an error when xleft and xright do not have different
  %signs
  %error prints the message in red and stops execution
  %sign(x): x > 0 -> sign(x) = 1, x < 0 -> sign(x) = -1, x = 0 ->
  %sign(x) = 0
  if sign(func(xleft)) == sign(func(xright))
    error('guesses must bracket zero')
  end
  
  while ~done
    iteration = iteration+1;
    fprintf('%2d:  %10.8f  %10.8f  %10.8f  %10.8f', ...
      iteration, xleft, func(xleft),xright, func(xright))

    xnew = (xleft + xright)/2;
    %Check if (approximate) error is below epsilon
    if abs(func(xnew)) < epsilon
      done = 1;
    elseif sign(func(xnew)) == sign(func(xleft))
      xleft = xnew;
    else
      xright = xnew;
    end

    fprintf('  %10.8f  %10.8f\n',  xnew, func(xnew))
  end

end

