The examples of ApproxFun.jl examples solutions are fascinating and powerful, and I’d love to figure out the best performance for solving the homogeneous screened Poisson equation but in cylindrical coordinates (with circular symmetry) and with a mix of Dirichlet and homogeneous Neumann boundaries.

Here with the particular cylindrical Laplacian and circular symmetry (\tfrac{\partial}{\partial \theta} = 0):

The ApproxFun approach below seems to solve this problem, albeit a bit slowly. It is slower when setting constant Dirichlet boundaries than when using the \exp(r)\cos(z) expressions. I’ve reduced the data using `resizedata!`

and reduced the tolerance with `tolerance=1E-5`

.

```
λ = 10
a = 0.5
k_1 = 1.2
k_2 = 0.2
k_3 = 1.0
dr = a .. 1.0
dz = 0.0 .. 1.0
d = dr × dz
r, z = Fun(∂(d))
r, z = components(r), components(z)
g = [-exp(r[1])*sin(0); exp(1)*cos(z[2]);
exp(r[3])*cos(1); exp(a)*cos(z[4])]
# g2 = [0; k_3; k_1; k_2]
Δ = Laplacian(d)
Dr = Derivative(d, [1, 0])
B = [I⊗lneumann(dz); rdirichlet(dr)⊗I; I⊗rdirichlet(dz); ldirichlet(dr)⊗I]
invfun = Fun((r, z) -> 1 / r, d)
L = [B; invfun * Dr + Δ - λ^2 * I]
Q = qr(L)
ApproxFunBase.resizedata!(Q,:,10);
u = \(Q, [g; 0]; tolerance=1E-5)
```

What else can I do to speed up the solution?

PS: Thank you for making these tools in the first place.