I identified several issues with the coef_func:
Basic setup:
using DiffEqOperators
nknots = 10
h = 1.0/nknots
x = 1.0:10.0
bc = DirichletBC(0.0, 11.0)
Δ = CenteredDifference(1, 2, h, nknots)
Δ * (bc * x)
This returns a vector of 10, as expected.
If I pass any callable to the coef_func argument of CenteredDifference, the operator totally broke down:
Δ = CenteredDifference(1, 2, h, nknots, x-> h)
Δ * (bc * x)
This returns a vector of 0.
Of course, the coefficient for centereddifference function doesn't really matter, but other issues arise for UpwindDifference, mainly related with the treatment of boundaries:
Δ = UpwindDifference(1, 1, h, nknots, x-> h * x)
Δ * (bc * x)
This returns [1,2,1,2,3,4,5,6,-9,-10]. It seems that the coefficients are mismatched: They are only applied from the third entry. Also, there seem to be a problem with the right boundary.
If we concretize Δ, the mismatching is fixed:
which returns [1,2,3,4,5,6,7,8,-9,-10], though there are still problems with the right boundary.