Interpolating with 3 Points (2nd Order)

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

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

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

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

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

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

>    simplify( answer );

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

Interpolating with 5 Points (4th Order)

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

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






>    derivative := diff( interpolating, x );

derivative := 1/6*1/h^4*(f(x0+2*h)-4*f(x0+h)+6*f(x0)+f(x0-2*h)-4*f(x0-h))*x^3+1/8*1/h^4*(4*h*f(x0-h)+2*f(x0+2*h)*h-4*f(x0+h)*h-2*f(x0-2*h)*h+16*f(x0+h)*x0+16*f(x0-h)*x0-4*f(x0-2*h)*x0-24*f(x0)*x0-4*f(x...




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

answer := 1/6*1/h^4*(f(x0+2*h)-4*f(x0+h)+6*f(x0)+f(x0-2*h)-4*f(x0-h))*x0^3+1/8*1/h^4*(4*h*f(x0-h)+2*f(x0+2*h)*h-4*f(x0+h)*h-2*f(x0-2*h)*h+16*f(x0+h)*x0+16*f(x0-h)*x0-4*f(x0-2*h)*x0-24*f(x0)*x0-4*f(x0+2...




>    simplify( answer );

1/12*1/h*(-f(x0+2*h)+8*f(x0+h)-8*f(x0-h)+f(x0-2*h))