I am using legacy FEniCS 2019.1.0 and I need to find the eigenvalues and eigenfunctions of a Matern covariance operator, which is defined using modified bessel functions. However, the following mwe, defining an expression using bessel functions, does not compile:

```
```from fenics import *

import numpy as np

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 20, 20)

V = FunctionSpace(mesh, ‘P’, 1)

u = TrialFunction(V)

v = TestFunction(V)

k = Expression(‘pow(sqrt(pow(x[0]-x_i,2)+pow(x[1]-y_i,2)),nu)*cyl_bessel_k(nu, sqrt(pow(x[0]-x_i,2)+pow(x[1]-y_i,2)))’, degree=1, nu=1, x_i=0, y_i=0)

c = np.zeros((mesh.num_vertices(), mesh.num_vertices()))

coords = mesh.coordinates()

dof2v = dof_to_vertex_map(V)

for i in range(len(coords)):

k.x_i = Vertex(mesh, dof2v[i]).midpoint()[0]

k.y_i = Vertex(mesh, dof2v[i]).midpoint()[1]

c[i, :] = assemble(k * v * dx)

vals, vecs = np.linalg.eig(c)

The error message is the following:

/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp: In member function ‘virtual void dolfin::dolfin_expression_4a2c2db2129efa596eb787cd10ed666d::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:

/tmp/tmp6vvyh3xk/dolfin_expression_4a2c2db2129efa596eb787cd10ed666d.cpp:64:69: error: ‘cyl_bessel_k’ was not declared in this scope

Apparently cyl_bessel_k is not declared even though it is listed in the documentation of cmath: Standard library header <cmath> - cppreference.com

I have also looked into this (old) solution but the code provided there does not compile either, giving an extremely long error message, mostly repeating error: ‘namespace’ definition is not allowed here https://fenicsproject.org/qa/3498/use-bessel-functions-in-jit-compiled-expressions/

Is there a different way of approaching this problem? Or am I just using the wrong alias for bessel functions?