Right, but I still need to construct the matrix first, even if I wrap it as an operator later. This is because the full matrix has some funky structure due to the underlying equation (see the second bullet point in the linked post). That summed V_{F0} term, combined with the fact that the operator spans all instances of the index i, means I need to explicitly construct the matrix operator for the PDE. I can do all the steps of the construction as long as I can get the 1D second-derivative operator on one axis of grid points, but doing that with FD methods for a non-uniform grid is where I’m getting hung up.