function [D, diffOp] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
    diffOp = scheme(g, order, schemeParams{:}, opSet);
    if ~isa(diffOp.D, 'function_handle')
        D = @(v) diffOp.D*v;
    else
        D = diffOp.D;
    end
end