function(f, a, b, tol, maxiters)
{

# bisection
# uses bisection algorithm to find root of func in interval [a,b]
# Burden and Faires, section 2.1

################################################################
# inputs
################################################################
# tol      -- maximum difference between subinterval endpoints to
#             consider root to have been found
# f        -- function for which root needs to be found
# a,b      -- interval endpoints, b > a
# maxiters -- maximum number of iterations
################################################################

# initial setup

if( f(a) * f(b) > 0)
	print("Function has same sign at both endpoints.")
else {
	absdiff <- b-a
	iters <- 1

	while ((absdiff > tol) & (iters <= maxiters))
        	{
        	absdiff <- b-a # note: absdiff is constructed to be positive
        	p <- a + absdiff / 2
        	if ((f(p) != 0) && (absdiff > tol)) {
			if( f(p) * f(a) < 0 )
				b <- p
			else
				a <- p
			iters <- iters + 1
		}
	}		
	if (absdiff > tol)   # didn't find solution in fewer than maxiters
		print("Maximum number of iterations exceeded.")

	list( a = a, b = b, p = p, errflag = as.numeric(absdiff > tol) )
}

}