I love the following example of solving a 1D Poisson equation with DiffEqOperators.

However , I struggle to do the analog 3D example. Can somebody indicate me the right steps?

using DiffEqOperators, Plots

#Poisson Eq in 1D

#`Δu = f`

with boundary conditions `u(0) = 0`

and `u(1) = 0`

Δx = 0.1

x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!

N = length(x)

f = sin.(2π*x)

Δ = CenteredDifference(2, 2, Δx, N)

bc = Dirichlet0BC(Float64)

u = (Δ*bc) \ f

u_analytic(x) = -sin(2π*x)/4(π^2)

plot(x, u)

plot!(x, u_analytic.(x))

In 3D, I could define Δ as:

Δxx = CenteredDifference{1}(2,2,Δx,N)

Δyy = CenteredDifference{2}(2,2,Δx,N)

Δzz = CenteredDifference{3}(2,2,Δx,N)

Δ = Δxx + Δyy + Δzz

But then, what to do ?

Thanks in advance for any kind of answer …