There used to be a nice FDFormula
at WRI site, but it is gone now. But I used it before. Here is the result.
I'll show some examples then the code at end
getFormula[1, {-1, 0, 1}, "centered"]

The first argument to getFormula
is derivative order. So 1
for first order, 2
for second order. The second argument is list of points to generate the difference approximation on. The last argument is the type you want. Either centered, forward or backward.
The function returns the difference formula and also the error in the approximation (the big O).
Here are more examples
getFormula[1, {-1, 0, 1}, "forward"]

getFormula[1, {-1, 0, 1}, "backward"]

Second order
getFormula[2, {-1, 0, 1}, "centered"]

More points, gives better approximations
getFormula[2, {-2, -1, 0, 1, 2}, "centered"]

getFormula[2, {-1, 0, 1}, "backward"]

4th order. Need to supply more grid points in this case, otherwise will get an error.
getFormula[4, {-2, -1, 0, 1, 2}, "centered"]

Code
(*FDFormula from
http://reference.wolfram.com/mathematica/tutorial/NDSolvePDE.html*)
FDFormula[(m_Integer)?Positive, (n_Integer)?Positive, (s_Integer)?
NonNegative] :=
Module[{do, F}, F = Table[f[Subscript[x, i + k]], {k, -s, n - s}];
W =
PadRight[
CoefficientList[Normal[Series[x^s*Log[x]^m, {x, 1, n}]/h^m], x],
Length[F], 0];
Wfact = 1/PolynomialGCD @@ W; W = Simplify[W*Wfact];
taylor[h_] =
Normal[Series[f[Subscript[x, i] + h], {h, 0, n + 2}]];
error = Drop[CoefficientList[
Expand[Table[taylor[h*k], {k, -s, n - s}] . W/Wfact], h], 1];
do = Position[error, e_ /; e != 0][[1, 1]]; error = error[[do]];
error = error /. (f_)[Subscript[x, i]] -> f; error = h^do*error;
{Derivative[m][f][Subscript[x, i]] \[TildeEqual] F . W/Wfact,
error}]
This uses the above function
getFormula[order_, gridPoints_, type_String] := Module[{s},
s = Which[type == "centered", (Length[gridPoints] - 1)/2,
type == "forward", 0,
True, Length[gridPoints] - 1];
Print[s];
FDFormula[order, Length[gridPoints] - 1, s]
]
I've used this in past to make a detailed Manipulate. But I never send it to Wolfram demo site.

This Demonstration illustrates the effect of numerical errors on the
approximation of derivatives when using the finite-difference scheme
with different step sizes and different orders of accuracy. You can
select to approximate up to the fourth derivative, the desired local
truncation accuracy order O(h^n), and the finite difference scheme to
use (centered, forward, or backward).
Also @xzczd has a finite difference formula generator function on this site. I do not have the link right now. That might also be something to look at. I've seen him use it to answer many questions.
ps. if you want to download the full Manipulate shown above, you can go to this page and search for "difference" and you'll find it near the top of the page there. One day I might submit to WRI demo site when I clean it a little more.
D[\[Phi]1, x, x]
should beD[\[Phi]1[x, y], x, x]
. 2.Replace[F, replacements]
should beReplace[F, replacements, 1]
, check the document ofReplace
for more info. Alternatively, useReplaceAll
(/.
) instead. 3.\[Delta]x
should be\[Delta]x^2
. $\endgroup$D[{ϕ1[x, y], ϕ2[x, y], ϕ3[x, y]}, {x, 2}] /. Derivative[ords__][f_][args__] :> Fold[(DifferenceQuotient[f[args], {#[[2]], #[[1]] - 1, h}] - DifferenceQuotient[f[args], {#[[2]], #[[1]] - 1, -h}])/h &, Transpose[{{ords}, {args}}]] // Simplify
. $\endgroup$