SkeelKeiper Improved Bisection (p. 79)

Printout[output_,a_,b_,c_,fc_,slope_]:=
Module[{},
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"(Subscript[a, k])","(Subscript[c, k])","(Subscript[b, k])","f[(Subscript[c, k])]","slope"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",((b-a)/2)]; 
Print["f[c] = ",NumberForm[fc,16] ]; 
Print["f'[c] = ",NumberForm[slope,16] ]; 
]

Options[Bisection]={maxits->100,epsilon->0} ;
Bisection[f_,a0_,b0_,OptionsPattern[]] :=
Module[{i=1,its=OptionValue[maxits],eps=OptionValue[epsilon]},
If[Abs[f[a0]]≤eps, Print["Root at endpoint a."];Return[Transpose[{{a0,a0,b0,0}}]]];
If[Abs[f[b0]]≤eps, Print["Root at endpoint b."];Return[Transpose[{{a0,b0,b0,0}}]]];
a=N[a0]; (* Slight change -- use a local variable (so as not to alter passed arguments *)
b=N[b0];
fa = f[a];  
fb = f[b];  
If[fa*fb>=0,Print["Error: function value doesn't change sign."];Return[Transpose[{{a0,a0,b0,0}}]]];
width=b-a;
slope=(fb-fa)/width;
c = a+0.5*width; 
fc=f[c]; (* Slight change -- moved first evaluation of fc above the while, immediately after c is defined.*)
output={{a,c,b,fc,slope}}; (* Let's collect the output, since we've gone to the trouble of computing it. *)
If[Abs[fc]≤eps, Print["Root at c."];Printout[output,a,b,c,fc,slope];Return[Transpose[{{a,c,b,fc}}]]];While[a < c && c<b  && i<its,
If[ Sign[fa] ==  Sign[fc], 
a = c;fa=fc, 
b = c;fb=fc
];
width=b-a;
slope =(fb-fa)/width;
i=i+1;
c = a+0.5*width;
fc=f[c];
If[Abs[fc]≤eps, Print["Root at c."];Printout[output,a,b,c,fc,slope];Return[Transpose[{{a,c,b,fc}}]]];
output=Append[output,{a,c,b,fc,slope}]; (* once again, collect output. *)
]; 
(* Now let's augment the naive function to include printing some nice formatted output: *)
Printout[output,a,b,c,fc,slope];
Return[Transpose[output]]; (* Now I'll actually return the output, so that it can be plotted, etc. *)
]

We first illustrate that the naive method works (to some extent!):

f[x_]:=Cos[x]-x
results=Bisection[f,0,Pi/2,maxits->10];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

What if our inputs don’t have function values of opposite sign?

results=Bisection[f,Pi/3,Pi/4,maxits->5];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Now show that it can be used to find quotients, using only multiplication, addition, and subtraction:

f[x_]:=3x-1
results=Bisection[f,.1,1];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Our interval enpoints are roots, and now we catch them:

f[x_]:=Sin[x]
results=Bisection[f,0,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]
results=Bisection[f,0.1,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Our first step is onto a root, but we miss it -- because of machine error! Our algorithm converts the answers to decimal, which has the effect of ruining a perfectly good root. Notice how long it takes us to get even somewhat close to the initial value of f[c] (our final result isn’t as good as our first calculated approximation).

f[x_]:=Cos[x]
results=Bisection[f,0,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

We can catch it by reducing our demand on f[c]: we give an epsilon of the appropriate size that we desire:

f[x_]:=Cos[x]
results=Bisection[f,0,Pi,epsilon->10^-6];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

What if we try to find a root of a discontinuous function?

f[x_]:=.5-Floor[x]
results=Bisection[f,-.1,Pi,epsilon->10^-6];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]