Finite difference operators for PDE solving in Julia (DiffEqOperators alternative)

Question

In the past I was using DiffEqOperators.jl to solve PDEs on a nonuniform grid, but now that this package has been archived I am struggling to find an alternative which provides the same functionality.

I have a Fokker-Planck-like equation of the form

du(x,t)/dt = d[B(x)u(x,t)]^2/dx^2 + C*delta(x-z)

which is made somewhat complicated by the delta function inside the boundaries.

I was using a method-of-lines approach as outlined by Chris Rackauckas in this overview https://nextjournal.com/sosiris-de/pde-2018. Basically, I would take a (nonuniform) grid dx_i (a 1D-vector of length L) upon which I could directly obtain the differential operator of desired order Δ = CenteredDifference(2, order, dx_i, L). With this I could calculate Δ*B*u_i for all points, to obtain the differential du_i .= Δ*B*u_i + C*deltaZ_i at every grid point (where deltaZ_i is some representation of the delta function at $z$). This I could then evolve with OrdinaryDiffEq.

If I understand correctly, the MethodOfLines.jl package is now recommended for this type of discretization and subsequent MOL evolution, however I've found it neither as performant nor as accurate as my previous "manual" implementation (though perhaps that might just be me misusing it), and would prefer to implement my own discretization.

My question is thus if there still is a package which implements finite difference operators -- on non-uniform grids and to (somewhat) arbitrary order -- that can be applied to a function for which I only have the values at my selected grid points dx_i.

I've looked at both FiniteDifferences.jl and FiniteDiff.jl, but these appear to require functional forms f(x), for which I would need to do an additional interpolation step.

For example, in FiniteDiff.jl to calculate the derivative one calls

FiniteDiff.finite_difference_derivative(
    f,
    x          :: AbstractArray{<:Number},
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x),      # return type of f
    fx         :: Union{Nothing,AbstractArray{<:Number}} = nothing,
    epsilon    :: Union{Nothing,AbstractArray{<:Real}} = nothing;
    [epsilon_factor])

(or potentially the in-place or caching form), which calculates the derivative on f, which must be a function; whereas my solution u_i only has values at gridpoints i. The same goes for the FiniteDifferences.jl package. Maybe I am missing something here?