diff --git a/demo/lbb-stability/lbb-stability.py b/demo/lbb-stability/lbb-stability.py
index af316ff..2942dbe 100644
--- a/demo/lbb-stability/lbb-stability.py
+++ b/demo/lbb-stability/lbb-stability.py
@@ -2,7 +2,9 @@
 import argparse
 import pathlib
 import firedrake
-from firedrake import assemble, inner, sym, grad, div, dx
+from firedrake import (
+    assemble, inner, sym, grad, div, dx, FiniteElement, TensorProductElement
+)
 import matplotlib.pyplot as plt
 import numpy as np
 from petsc4py import PETSc
@@ -17,38 +19,38 @@ def compute_inf_sup_constant(spaces, b, inner_products):
     u, p = firedrake.TrialFunctions(Z)
     v, q = firedrake.TestFunctions(Z)
 
-    A = assemble(-(m(u, v) + b(v, p) + b(u, q)), mat_type="aij").M.handle
-    M = assemble(n(p, q), mat_type="aij").M.handle
-
-    opts = PETSc.Options()
-    opts.setValue("eps_gen_hermitian", None)
-    opts.setValue("eps_target_real", None)
-    opts.setValue("eps_smallest_magnitude", None)
-    opts.setValue("st_type", "sinvert")
-    opts.setValue("st_ksp_type", "preonly")
-    opts.setValue("st_pc_type", "lu")
-    opts.setValue("st_pc_factor_mat_solver_type", "mumps")
-    opts.setValue("eps_tol", 1e-8)
-
-    num_values = 1
-    eigensolver = SLEPc.EPS().create(comm=firedrake.COMM_WORLD)
-    eigensolver.setDimensions(num_values)
-    eigensolver.setOperators(A, M)
-    eigensolver.setFromOptions()
-    eigensolver.solve()
-
-    Vr, Vi = A.getVecs()
-    λ = eigensolver.getEigenpair(0, Vr, Vi)
-    return λ
-
-
-def run(*, dimension, num_cells, element, num_layers=None, vdegree=None):
+    A = -(m(u, v) + b(v, p) + b(u, q))
+    M = n(p, q)
+
+    problem = firedrake.LinearEigenproblem(A=A, M=M, restrict=False)
+    params = {
+        "n_evals": 1,
+        "solver_parameters": {
+            "eps_tol": 1e-8,
+            "eps_gen_hermitian": None,
+            "eps_target_real": None,
+            "eps_smallest_magnitude": None,
+            "st_type": "sinvert",
+            "st_ksp_type": "preonly",
+            "st_pc_type": "lu",
+            "st_pc_factor_mat_solver_type": "mumps",
+        },
+    }
+    solver = firedrake.LinearEigensolver(problem, **params)
+    num_converged = solver.solve()
+    assert num_converged >= 1
+    return solver.eigenvalue(0)
+
+
+def run(
+    *, dimension, num_cells, element, num_layers=None, vdegree=None, vdegree_delta=None
+):
     mesh = firedrake.UnitSquareMesh(num_cells, num_cells, diagonal="crossed")
     h = mesh.cell_sizes.dat.data_ro[:].min()
 
-    cg_1 = firedrake.FiniteElement("CG", "triangle", 1)
-    cg_2 = firedrake.FiniteElement("CG", "triangle", 2)
-    b_3 = firedrake.FiniteElement("Bubble", "triangle", 3)
+    cg_1 = FiniteElement("CG", "triangle", 1)
+    cg_2 = FiniteElement("CG", "triangle", 2)
+    b_3 = FiniteElement("Bubble", "triangle", 3)
     if element == "taylor-hood":
         velocity_element = cg_2
         pressure_element = cg_1
@@ -57,13 +59,13 @@ def run(*, dimension, num_cells, element, num_layers=None, vdegree=None):
         pressure_element = cg_1
     elif element == "crouzeix-raviart":
         velocity_element = cg_2 + b_3
-        pressure_element = firedrake.FiniteElement("DG", "triangle", 1)
+        pressure_element = FiniteElement("DG", "triangle", 1)
 
     if dimension == 3:
-        cg_z = firedrake.FiniteElement("CG", "interval", vdegree)
-        dg_z = firedrake.FiniteElement("DG", "interval", vdegree - 1)
-        velocity_element = firedrake.TensorProductElement(velocity_element, cg_z)
-        pressure_element = firedrake.TensorProductElement(pressure_element, dg_z)
+        cg_z = FiniteElement("CG", "interval", vdegree)
+        dg_z = FiniteElement("DG", "interval", vdegree - vdegree_delta)
+        velocity_element = TensorProductElement(velocity_element, cg_z)
+        pressure_element = TensorProductElement(pressure_element, dg_z)
         mesh = firedrake.ExtrudedMesh(mesh, num_layers)
 
     V = firedrake.VectorFunctionSpace(mesh, velocity_element)
@@ -95,6 +97,7 @@ def run(*, dimension, num_cells, element, num_layers=None, vdegree=None):
 parser.add_argument("--num-cells-step", type=int, default=1)
 parser.add_argument("--num-layers", type=int, default=1)
 parser.add_argument("--vdegree", type=int, default=1)
+parser.add_argument("--vdegree-delta", type=int, default=1)
 parser.add_argument("--output")
 args = parser.parse_args()
 
@@ -109,7 +112,8 @@ def run(*, dimension, num_cells, element, num_layers=None, vdegree=None):
         num_cells=num_cells,
         element=args.element,
         num_layers=args.num_layers,
-        vdegree=args.vdegree
+        vdegree=args.vdegree,
+        vdegree_delta=args.vdegree_delta,
     )
     stability_constants.append((h, λ))
     print(".", flush=True, end="")