The derivative operators are internally converted to sparse matrices anyway, so your approach is in principle good.
Beware that one needs to solve (I + Dxx*bc) \ (b-bxx) though. See https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf for a discussion of ghost nodes and boundary conditions, and take a look at how DiffEqOperators implements \.