Here is a method that might be more automated than the one
you indicated. In wxmaxima, positional derivatives
display as f[(1)], f[(2)], ..., where f([n]), is the n-th
derivative of f with respect to the first argument.
(%i1) load("pdiff")$
(%i2) eq : f(x+a) = taylor(f(x+a),a,0,4)$
(%i3) eqs : makelist(subst(a = k*h,eq),k,-2,2)$
(%i4) linsolve(eqs, makelist(diff(f(x),x,k),k,1,4));
Dependent equations eliminated: (3)
(%o4) [f[(1)](x)=(-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h),f[(2)](x)=-
(f(x+2*h)-16*f(x+h)-16
*f(x-h)+f(x-2*h)+30*f(x))/(12*h^2),f[(3)](x)=-(-f(x+2*h)+2*f(x+h)-2
*f(x-h)+f(x-2*h))/(2*h^3),f[(4)](x)=
(f(x+2*h)-4*f(x+h)-4*f(x-h)+f(x-2*h)+6*f(x))/h^4]
Look at the leading order of the error terms:
(%i7) map(lambda([e], lhs(e) = taylor(rhs(e),h,0,5)),%o4);
(%o7)
[f[(1)](x)=f[(1)](x)-(f[(5)](x)*h^4)/30+...,f[(2)](x)=f[(2)](x)-(f[(6)](x)*h^4)/90+...,f[(3)](x)=f[(3)](x)+(f[(5)](x)*h^2)/4+
(f[(7)](x)*h^4)/40+...,f[(4)](x)=f[(4)](x)+(f[(6)](x)*h^2)/6+(f[(8)](x)*h^4)/80+...]
(%i8)
A similar method works for FD approximations for functions of two or
more variables.
Barton