FDFormula[m_(*order of derivative*), n_(*number of grid intervals*),
s_(*set to n/2 to get centered differenrer*)] :=
Block[{}, F = Table[f[Subscript[x, i + k]], {k, -s, n - s}];
W = PadRight[
CoefficientList[Normal[Series[x^s Log[x]^m, {x, 1, n}]/h^m], x],
Length[F], 0];
Wfact = 1/Apply[PolynomialGCD, W];
W = Simplify[W Wfact];
taylor[h_] = Normal[Series[f[Subscript[x, i] + h], {h, 0, n + 2}]];
error =
Drop[CoefficientList[
Expand[(Table[taylor[h k], {k, -s, n - s}].W)/Wfact], h], 1];