Demo that Bisection doesn’t converge linearly

This just defines bisection:

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. *)
]

Here’s the demo:

This is a simple function, with root 4/5: starting from the interval [0,1], we will demonstrate non-linear convergence....

f[x_]:=x-4/5
results=Bisection[f,0,1];
cs=results[[2]]; (* the computed c values *)
n=Length[cs];

Now we ttake the absolute values of the errors, and look at the ratios e_{k+1}/e_k: they oscillate (and don’t tend to a fixed limit):

errors=Abs[cs-4/5]
ListPlot[errors[[2;;n]]/errors[[1;;n-1]]]