CDDF stands for Cendered Divided Difference Formula
> 
CDDF := proc(f, A, H, n, stencil)
description "f = target function,\nA = point where the derivatives are approximated,\nH = step,\nn = order of the derivative,\nstencil = list of points for the divided differenceCDDF\n";
local tay, p, T, sol, unknown, Unknown, expr:
tay := (s, m) > convert(
eval(
convert(
taylor(op(0, f)(op(1, f)), op(1, f)=A, m),
Diff
),
op(1, f)=A+s*H),
polynom
):
p := numelems(stencil):
T := add(alpha[i]*tay(i, p+1), i in stencil):
T := convert(%, diff):
if p > n+1 then
sol := solve([seq(coeff(T, h, i)=0, i in subsop(n+1=NULL, [$0..p]))], [seq(alpha[i], i in stencil)])[];
else
sol := solve([seq(coeff(T, H, i)=0, i in subsop(n+1=NULL, [$0..n]))], [seq(alpha[i], i in stencil)])[];
end if:
if `and`(is~(rhs~(sol)=~0)[]) then
WARNING("no solution found"):
return
else
unknown := `union`(indets~(rhs~(sol))[])[];
Unknown := simplify(solve(eval(T, sol) = Diff(op(0, f)(A), A$n), unknown));
sol := lhs~(sol) =~ eval(rhs~(sol), unknown=Unknown);
expr := normal(eval(add(alpha[i]*op(0, f)(A+i*H), i in stencil), sol));
end if:
return expr
end proc:

# f = target function,
# A = point where the derivatives are approximated,
# H =
# step,
# n = order of the derivative,
# stencil = list of points for the divided
# differenceCDDF
#
CDDF( f, A, H, n, stencil )


> 
# 2point approximation of diff(f(x), x)  x=a
CDDF(f(x), a, h, 1, [1, 1]);
convert(simplify(mtaylor(%, h=0, 4)), Diff);


(1) 
> 
# 3point approximation of diff(f(x), x$2)  x=a
CDDF(f(x), a, h, 2, [1, 0, 1]);
convert(simplify(mtaylor(%, h=0)), Diff);


(2) 
> 
# 5point pproximation of diff(f(x), x$2)  x=a
CDDF(f(x), a, h, 2, [$2..2]);
convert(simplify(mtaylor(%, h=0, 8)), Diff);


(3) 
> 
# 7point approximation of diff(f(x), x$2)  x=a
CDDF(f(x), a, h, 2, [$3..3]);
# simplify(taylor(%, h=0, 10));
convert(simplify(mtaylor(%, h=0, 10)), Diff);


(4) 
> 
# 4point staggered approximation of diff(f(x), x$3)  x=a
CDDF(f(x), a, h, 3, [seq(3/2..3/2, 1)]);
convert(simplify(mtaylor(%, h=0, 6)), Diff);


(5) 
> 
# 6point staggered approximation of diff(f(x), x$3)  x=a
CDDF(f(x), a, h, 3, [seq(5/2..5/2, 1)]);
# simplify(taylor(%, h=0, 8));
convert(simplify(mtaylor(%, h=0, 8)), Diff);


(6) 
> 
# 5point approximation of diff(f(x), x$4)  x=a
CDDF(f(x), a, h, 4, [$2..2]);
convert(simplify(mtaylor(%, h=0, 8)), Diff);


(7) 
> 
# 7point approximation of diff(f(x), x$4)  x=a
CDDF(f(x), a, h, 4, [$3..3]);
convert(simplify(mtaylor(%, h=0, 10)), Diff);


(8) 
A FEW 2D EXTENSIONS
> 
# diff(f(x, y), [x, y]) approximation over a (2 by 2)point stencil
stencil := [1, 1]:
# step 1: approximate diff(f(x, y), x) over stencil < stencil >
fx := CDDF(f(x), a, h, 1, stencil):
fx := eval(% , f=(u > f[u](y))):
ix := [indets(fx, function)[]]:
# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.
fxy := add(map(u > CDDF(u, b, k, 1, stencil)*coeff(fx, u), ix)):
# step 3: rewrite fxy in a more convenient form
[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );
convert(mtaylor(fxy, [h=0, k=0]), Diff)


(9) 
> 
# Approximation of diff(f(x, y), [x, x, y, y] a (3 by 3)point stencil
stencil := [1, 0, 1]:
# step 1: approximate diff(f(x, y), x) over stencil < stencil >
fx := CDDF(f(x), a, h, 2, stencil):
fx := eval(% , f=(u > f[u](y))):
ix := [indets(fx, function)[]]:
# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.
fxy := add(map(u > CDDF(u, b, k, 2, stencil)*coeff(fx, u), ix)):
# step 3: rewrite fxy in a more convenient form
[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );
convert(mtaylor(fxy, [h=0, k=0], 8), Diff)


(10) 
> 
# Approximation of diff(f(x, y), [x, x, y] a (3 by 2)point stencil
stencil_x := [1, 0, 1]:
stencil_y := [1, 1]:
# step 1: approximate diff(f(x, y), x) over stencil < stencil >
fx := CDDF(f(x), a, h, 2, stencil_x):
fx := eval(% , f=(u > f[u](y))):
ix := [indets(fx, function)[]]:
# step 2: approximate diff(g(y), y) over stencil < stencil > where
# g represents any function in fx.
fxy := add(map(u > CDDF(u, b, k, 1, stencil_y)*coeff(fx, u), ix)):
# step 3: rewrite fxy in a more convenient form
[seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]:
fxy := simplify( eval(fxy, %) );
convert(mtaylor(fxy, [h=0, k=0], 6), Diff)


(11) 
> 
# Approximation of the laplacian of f(x, y)
stencil := [1, 0, 1]:
# step 1: approximate diff(f(x, y), x) over stencil < stencil >
fx := CDDF(f(x), a, h, 2, stencil):
fy := CDDF(f(y), b, k, 2, stencil):
fxy := simplify( eval(fx, f=(u > f(u, b))) + eval(fy, f=(u > f(a, u))) );
convert(mtaylor(fxy, [h=0, k=0], 6), Diff)


(12) 
