I need to construct a replacement rule to replace (first) derivatives with centered finite differences. What I've got so far is
DtoFD := f_'[var_] -> (f[var + ϵ/2] - f[var - ϵ/2])/ϵ
which works on simple cases like
f'[x] /. DtoFD (* (-f[x - ϵ/2] + f[x + ϵ/2])/ϵ *)
Unfortunately I didn't consider the possibility of multivariate functions like
D[f[x, y], x] /. DtoFD (* f^(1,0)[x, y] *)
where it fails.
Could someone come up with a rule that works in the multivariate case? First differences would be enough, but I'll also take arbitrary order if you've got it!
EDIT: I should add my motivation is to automatically set up a Jacobian matrix for a system of ordinary difference / differential equations, where some of the right hand sides might be numerically defined "black box" functions. So this substitution rule serves as a backup in case some derivatives can't be found with