Numerical Differentiation using Mathematica

I am teaching myself Mathematica over summer, and coded up a simple module that computes numerical differentiation formulas automatically. Essentially, the module enables one to compute arbitrary differentiation rules like those explained here.

This module spits out the different formulae for a m-point approximation to the n-th derivative, and the leading error term.

Here is the (updated) module:

AppDeriv[m_, n_] := Module[{points, IP, nthDeriv, e1},
  points = Table[{x0 + i h, Subscript[f, i]}, {i, 0, m - 1}];
  IP = InterpolatingPolynomial[points, x];
  nthDeriv = D[IP, {x, n}];
  e1 = (1/Factorial[m])*
    D[Apply[Times, (x - Table[{x0 + i h}, {i, 0, m - 1}])], {x, n}];
  Formula =
   TableForm[
    Table[{Subsuperscript[f, i, n], nthDeriv /. x -> (x0 + i h),
       Superscript[f, m] e1 /. x -> (x0 + i h)}, {i, 0, m - 1}] //
     Simplify ]]

This seems to work alright for me (click to enlarge).