FiniteDifferenceFormula:=proc(m,k,j)
# computes for m+1 equidistant points -jh, ..., -h, 0, h, ...,(m-j)h
# a finite difference approximation of the kth derivative evaluated 
# at x=0
  local i, p, dp;
  p:=interp([seq(i*h,i=-j..(m-j))],[seq(f(i*h),i=-j..(m-j))],x);
  dp:=diff(p,x$k);
  simplify(eval(dp,x=0));
end :