diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml
index c0682518a4..b1cf6132b6 100644
--- a/.github/workflows/core.yml
+++ b/.github/workflows/core.yml
@@ -214,6 +214,8 @@ jobs:
             --extra-index-url https://download.pytorch.org/whl/cpu \
             "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]"
 
+          pip install -I git+https://github.com/firedrakeproject/ufl.git@pbrubeck/hypergeometric
+          pip install -I git+https://github.com/pbrubeck/loopy.git@pbrubeck/hypergeometric
           firedrake-clean
           pip list
 
diff --git a/pyop2/global_kernel.py b/pyop2/global_kernel.py
index 4dd404b1f3..c9c294aab3 100644
--- a/pyop2/global_kernel.py
+++ b/pyop2/global_kernel.py
@@ -396,7 +396,7 @@ def _cppargs(self):
     def _ldargs(self):
         ldargs = [f"-L{d}/lib" for d in get_petsc_dir()]
         ldargs.extend(f"-Wl,-rpath,{d}/lib" for d in get_petsc_dir())
-        ldargs.extend(["-lpetsc", "-lm"])
+        ldargs.extend(["-lpetsc", "-lm", "-lgsl", "-lgslcblas"])
         ldargs.extend(self.local_kernel.ldargs)
         return tuple(ldargs)
 
diff --git a/tests/firedrake/regression/test_hypergeometric_function.py b/tests/firedrake/regression/test_hypergeometric_function.py
new file mode 100644
index 0000000000..20c8a7e205
--- /dev/null
+++ b/tests/firedrake/regression/test_hypergeometric_function.py
@@ -0,0 +1,23 @@
+from firedrake import *
+from scipy.special import hyp2f1 as sp_hyp2f1
+import numpy as np
+import pytest
+
+
+def test_hypergeometric_function():
+    mesh = UnitDiskMesh(3)
+    V = FunctionSpace(mesh, "CG", 1)
+    u = Function(V)
+    w = Function(V)
+    x, y = SpatialCoordinate(mesh)
+
+    expressions = [((1/3, 1/2, 1), Constant(0.9999) * x),
+                   ((1/2, 1/2, 1), Constant(0.9999) * sqrt(x**2+y**2)),
+                   ((-1/2, 1/2, 1), Constant(0.4999) * sqrt(x*x + y*y))]
+
+    for (a, b, c), expr in expressions:
+        u.interpolate(hyp2f1(a, b, c, expr))
+        result = u.dat.data
+        w.interpolate(expr)
+        expect = sp_hyp2f1(a, b, c, w.dat.data)
+        assert np.allclose(result, expect)
diff --git a/tsfc/loopy.py b/tsfc/loopy.py
index 6826f0b672..88c7e2403e 100644
--- a/tsfc/loopy.py
+++ b/tsfc/loopy.py
@@ -480,7 +480,11 @@ def _expression_mathfunction(expr, ctx):
                 else:
                     return p.Variable(f"{name}n")(nu_, arg_)
     else:
-        if expr.name == "ln":
+        if expr.name == "hyp2f1":
+            assert isinstance(ctx.target, lp.target.c.CWithGNULibcTarget)
+            # Generate right functions calls to gsl hypergeometric function
+            name = "hyperg_2F1"
+        elif expr.name == "ln":
             name = "log"
         else:
             name = expr.name
diff --git a/tsfc/ufl2gem.py b/tsfc/ufl2gem.py
index e283fb1414..19e4967d17 100644
--- a/tsfc/ufl2gem.py
+++ b/tsfc/ufl2gem.py
@@ -81,6 +81,9 @@ def math_function(self, o, expr):
     def atan2(self, o, y, x):
         return MathFunction("atan2", y, x)
 
+    def hypergeometric2_f1(self, o, a, b, c, arg):
+        return MathFunction(o._name, a, b, c, arg)
+
     def bessel_i(self, o, nu, arg):
         return MathFunction(o._name, nu, arg)