Interpolating with 2 Points

>    interpolating := CurveFitting[PolynomialInterpolation]( [x0 - h, x0], [f(x0 - h), f(x0)], x );

interpolating := (f(x0)-f(x0-h))/h*x+(f(x0-h)*x0-f(x0)*x0+f(x0)*h)/h

>    derivative := diff( interpolating, x ); # differentiate w.r.t. x

derivative := (f(x0)-f(x0-h))/h

>    answer := eval( derivative, x = x0 ); # evaluate at x0

answer := (f(x0)-f(x0-h))/h

>    simplify( answer );

(f(x0)-f(x0-h))/h

Interpolating with 3 Points

>    interpolating := CurveFitting[PolynomialInterpolation]( [x0 - 2*h, x0 - h, x0], [f(x0 - 2*h), f(x0 - h), f(x0)], x );

interpolating := 1/2*(f(x0)-2*f(x0-h)+f(x0-2*h))/h^2*x^2+1/2*(f(x0-2*h)*h+4*f(x0-h)*x0-4*f(x0-h)*h-2*f(x0)*x0+3*f(x0)*h-2*f(x0-2*h)*x0)/h^2*x+1/2*(-3*f(x0)*x0*h+f(x0-2*h)*x0^2-f(x0-2*h)*x0*h+f(x0)*x0^2...

>    derivative := diff( interpolating, x );

derivative := (f(x0)-2*f(x0-h)+f(x0-2*h))/h^2*x+1/2*(f(x0-2*h)*h+4*f(x0-h)*x0-4*f(x0-h)*h-2*f(x0)*x0+3*f(x0)*h-2*f(x0-2*h)*x0)/h^2

>    answer := eval( derivative, x = x0 );

answer := (f(x0)-2*f(x0-h)+f(x0-2*h))/h^2*x0+1/2*(f(x0-2*h)*h+4*f(x0-h)*x0-4*f(x0-h)*h-2*f(x0)*x0+3*f(x0)*h-2*f(x0-2*h)*x0)/h^2

>    simplify( answer );

1/2*1/h*(f(x0-2*h)-4*f(x0-h)+3*f(x0))