"""
c = mkfdstencil(x, xbar, k)
Compute the coefficients `c` in a finite difference approximation of a function
defined at the grid points `x`, evaluated at `xbar`, of order `k`.
"""
function mkfdstencil(x, xbar, k)
n = length(x)
A = @. (x[:]' - xbar) ^ (0:n-1) / factorial(0:n-1)
b = (1:n) .== k+1
c = A \ b
end
mkfdstencil
mkfdstencil([-1 0 1], 0, 2) # Centered 2nd derivative
3-element Vector{Float64}: 1.0 -2.0 1.0
mkfdstencil([0 1 2], 0, 1) # One-sided (right) 1st derivative
3-element Vector{Float64}: -1.5 2.0 -0.5
mkfdstencil(-2:2, 0//1, 2) # 4th order 5-point stencil, rational
5-element Vector{Rational{Int64}}: -1//12 4//3 -5//2 4//3 -1//12