I just started to play around with DiffEqOperators. I’d like to extend the example at

https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/examples/poisson.jl to a 2D case.

Here is the adapted code:

```
using DiffEqOperators, Plots
f = 5.0
a = -1.0
b = 2.0
nknots = 10
h = 1.0/(nknots+1)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
bc = DirichletBC(a, b)
u = (Δ*bc) \ fill(f, (nknots, nknots))
knots = range(h, step=h, length=nknots)
plot(knots, u[nknots÷2, 1:end])
plot!(knots, u[1:end, nknots÷2])
```

It works somehow - but so to say in 1D

So, my question: How do I define the BC for 2D (or more)? Or maybe there are more things to be changed in the script?