Linear Systems for Difference Formulas (Illustrating the idea of Richardson extrapolation, per p. 260, etc.) -- based on Taylor Series

A variety of difference formulas for approximating derivatives using the idea of the Taylor series:

Forward (or backward) difference formulas:

Clear[h]

This first one is the forward difference formula, p. 256:

TaylorCoefs=Transpose[
{
{1,0},
{1,h}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1}]

This one is the backward difference formula:

TaylorCoefs=Transpose[
{
{1,-h},
{1,0}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1}]

This is the centered difference formula, p. 257:

TaylorCoefs=Transpose[
{
{1,-h,h^2/2},
{1,0,0},
{1,h,h^2/2}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1, 0}]
Ans.{-h^3/3!,0,h^3/3!} (* order of accuracy *)

This one is a three-point (forward version of the backward three-point method at the top of 263):

TaylorCoefs=Transpose[
{
{1,0,0},
{1,h,h^2/2},
{1,2*h,(2*h)^2/2}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1, 0}]
Ans.{0,h^3/6,8h^3/6} (* order of accuracy *)

Here is a five-point centered method:

TaylorCoefs=Transpose[
{
{1,-2*h,(2*h)^2/2,-(2*h)^3/6,(2*h)^4/24},
{1,-h,h^2/2,-h^3/6,h^4/24},
{1,0,0,0,0},
{1,h,h^2/2,h^3/6,h^4/24},
{1,2*h,(2*h)^2/2,(2*h)^3/6,(2*h)^4/24}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1, 0,0,0}]
Ans.{-(2*h)^5/Factorial[5],-h^5/Factorial[5],0,h^5/Factorial[5],(2*h)^5/Factorial[5]}

Here is a five-point forward method:

TaylorCoefs=Transpose[
{
{1,0,0,0,0},
{1,h,h^2/2,h^3/6,h^4/24},
{1,2*h,(2*h)^2/2,(2*h)^3/6,(2*h)^4/24},
{1,(3*h),(3*h)^2/2,(3*h)^3/6,(3*h)^4/24},
{1,(4*h),(4*h)^2/2,(4*h)^3/6,(4*h)^4/24}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,1, 0,0,0}]
Ans.{0,h^5/Factorial[5],(2*h)^5/Factorial[5],(3*h)^5/Factorial[5],(4*h)^5/Factorial[5]}

We can create difference formulas for higher derivatives using the same idea of Taylor series:

Here is the centered difference approximation for the second derivative:

TaylorCoefs=Transpose[
{
{1,-h,h^2/2},
{1,0,0},
{1,h,h^2/2}
}];
SchemeSolve=LinearSolve[TaylorCoefs];
Ans=SchemeSolve[{0,0, 1}]
Ans.{h^4/Factorial[4],0,h^4/Factorial[4]}