Thank you @bhaveshshrimali. Well, I’m a bit overwhelmed by the requirements in order to install Dolphinx. I decided to try the other way around instead. Indeed, I’m dealing with the Helmholtz problem

\nabla ^2 p + \kappa^2 p = 0

for which the weak form is

\int_{\Omega}(\nabla p \cdot \nabla \overline{q} - \kappa^2 p \overline{q})d \Omega = \int_{\Gamma}(g \overline{q})ds

By looking through the provided link I wrote the following, where

e^{-i \kappa x_1} = cos(\kappa x_1) - i sin(\kappa x_1)

```
from fenics import *
import matplotlib.pyplot as plt
nx = 200
ny = 200
k = 15
#mesh and function space
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1*P1)
tol = 1E-14
def LEFT(x, on_boundary):
return on_boundary and abs(x[0] + 1.0) < tol
bound_x_re = Expression(("cos(k*x[0])"), degree=1, k=k)
bound_x_im = Expression(("-sin(k*x[0])"), degree=1, k=k)
#boundary conditions
bcs = [DirichletBC(V.sub(0), bound_x_re, LEFT),
DirichletBC(V.sub(1), bound_x_im, LEFT)]
#variational problem
(p_re, p_im) = TrialFunction(V)
(q_re, q_im) = TestFunction(V)
g_re = Constant(0.0)
g_im = Constant(0.0)
a = inner(grad(p_re),grad(q_re))*dx + inner(grad(p_im),grad(q_im))*dx - k*k*p_re*q_re*dx - k*k*p_im*q_im*dx
L = dot(g_re, q_re)*ds + dot(g_im, q_im)*ds
#compute solution
p = Function(V)
solve(a == L, p, bcs)
p_re, p_im = p.split()
plot(p_re)
plt.show()
```

which gives

Could you tell me if there’s any issue in my implementation?