This first set of inputs computes a table of center difference approximations with the value of h cut in half at each step.

Clear[f, x0, d, n1, n2];
f[x_] = Sin[x];
x0 = 1. ;
d[x_, h_] := (f[x0 + h] - f[x0 - h])/2/
        h
h = 12.8;
n1 = 
        Table[h = h/2; d[x0, h], {n, 1, 16}];
n1 // TableForm

Next, we improve the results from above using the first-order Richardson extrapolation.

n2 = 
        Table[n1[[n]] + (n1[[n]] - n1[[n - 1]])/3, {n, 2, 
            Length[n1]}];
n2 // TableForm

We can check if we have hit an appropriately small value of h by looking at the ratio ((N(h/2) - 
                N(h))/(N(h/4) - N(h)/
                  2)) .   Since the first-order Richardson extrapolation for center-differencing is O((h^2)), we expect the ratio to approach (2^2) = 4 as we cut h in half.  Here is a table of these ratios for the results from above:

Table[(n1[[n]] - n1[[n - 1]])/(n1[[n + 1]] - 
              n1[[n]]), {n, 2, Length[n1] - 1}] // TableForm

For comparison, here is a table of differences between the approximations in n2 and Mathematica's value for the derivative:

n2 - f'[x0] // TableForm