diff --git a/nutils/SI.py b/nutils/SI.py
index d0e0a627e..a86104001 100644
--- a/nutils/SI.py
+++ b/nutils/SI.py
@@ -391,6 +391,13 @@ def __stack_like(op, *args, **kwargs):
             raise TypeError(f'incompatible arguments for {op.__name__}: ' + ', '.join(dim.__name__ for dim in dims))
         return dims[0].wrap(op(arg0, *args[1:], **kwargs))
 
+    @register('prod', 'product')
+    def __prod(op, a, axis=None, *args, **kwargs):
+        (dim, arg), = Quantity.__unpack(a)
+        axes = range(arg.ndim) if axis is None else axis if isinstance(axis, tuple) else (axis,)
+        n = functools.reduce(operator.mul, [arg.shape[axis] for axis in axes], 1)
+        return (dim**n).wrap(op(arg, axis, *args, **kwargs))
+
     @register('curvature')
     def __evaluate(op, *args, **kwargs):
         (dim0, arg0), = Quantity.__unpack(args[0])
diff --git a/nutils/export.py b/nutils/export.py
index 8ec7034a0..7471a4a32 100644
--- a/nutils/export.py
+++ b/nutils/export.py
@@ -26,17 +26,13 @@ def mplfigure(name, /, **kwargs):
         :class:`matplotlib.figure.Figure` object.
     '''
 
-    import matplotlib.figure
-    import matplotlib.backends.backend_agg
-    fig = matplotlib.figure.Figure(**kwargs)
+    with warnings.ignore(DeprecationWarning):
+        from matplotlib.figure import Figure
+    fig = Figure(**kwargs)
     with log.userfile(name, 'wb') as f:
         yield fig
-        if f:
-            matplotlib.backends.backend_agg.FigureCanvas(fig)  # sets reference via fig.set_canvas
-            try:
-                fig.savefig(f, format=os.path.splitext(name)[1][1:])
-            finally:
-                fig.set_canvas(None)  # break circular reference
+        fig.savefig(f, format=os.path.splitext(name)[1][1:])
+        fig.clear()
 
 
 def plotlines_(ax, xy, lines, **kwargs):
diff --git a/nutils/function.py b/nutils/function.py
index 05e8aa073..d1937cd30 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -353,7 +353,7 @@ def prod(self, __axis: int) -> 'Array':
         :class:`Array`
         '''
 
-        return numpy.product(self, __axis)
+        return numpy.prod(self, __axis)
 
     def dot(self, __other: IntoArray, axes: Optional[Union[int, Sequence[int]]] = None) -> 'Array':
         '''Return the inner product of the arguments over the given axes, elementwise over the remanining axes.
@@ -3112,8 +3112,8 @@ def sum(arg: IntoArray, axis: Optional[Union[int, Sequence[int]]] = None) -> Arr
             summed = _Wrapper(evaluable.Sum, summed, shape=summed.shape[:-1], dtype=summed.dtype)
         return summed
 
-    @implements(numpy.product)
-    def product(arg: IntoArray, axis: int) -> Array:
+    @implements(numpy.prod)
+    def prod(arg: IntoArray, axis: int) -> Array:
         arg = Array.cast(arg)
         if arg.dtype == bool:
             arg = arg.astype(int)
@@ -3123,6 +3123,9 @@ def product(arg: IntoArray, axis: int) -> Array:
             multiplied = _Wrapper(evaluable.Product, multiplied, shape=multiplied.shape[:-1], dtype=multiplied.dtype)
         return multiplied
 
+    if hasattr(numpy, 'product'): # numpy < 2.0
+        implements(numpy.product)(prod)
+
     @implements(numpy.conjugate)
     def conjugate(arg: IntoArray) -> Array:
         return _Wrapper.broadcasted_arrays(evaluable.conjugate, arg)
diff --git a/nutils/matrix/_scipy.py b/nutils/matrix/_scipy.py
index 3c4270b05..b043f0a46 100644
--- a/nutils/matrix/_scipy.py
+++ b/nutils/matrix/_scipy.py
@@ -82,7 +82,8 @@ def mycallback(arg):
                 if callback:
                     callback(res)
                 reformat(100 * numpy.log10(max(atol, res)) / numpy.log10(atol))
-            lhs, status = solverfun(self.core, rhs, M=precon, tol=0., atol=atol, callback=mycallback, **solverargs)
+            solverargs['rtol' if scipy.version.version >= '1.12' else 'tol'] = 0.
+            lhs, status = solverfun(self.core, rhs, M=precon, atol=atol, callback=mycallback, **solverargs)
         if status != 0:
             raise Exception('status {}'.format(status))
         return lhs
diff --git a/nutils/testing.py b/nutils/testing.py
index c813f8068..451a97261 100644
--- a/nutils/testing.py
+++ b/nutils/testing.py
@@ -184,6 +184,7 @@ def setUp(self):
         self.enter_context(treelog.set(treelog.LoggingLog('nutils')))
         self.enter_context(_builtin_warnings.catch_warnings())
         _builtin_warnings.simplefilter('error', warnings.NutilsWarning)
+        _builtin_warnings.simplefilter('error', DeprecationWarning)
 
     def assertAllEqual(self, actual, desired):
         actual = numpy.asarray(actual)
diff --git a/nutils/warnings.py b/nutils/warnings.py
index d774ea8f3..0ab709911 100644
--- a/nutils/warnings.py
+++ b/nutils/warnings.py
@@ -31,4 +31,11 @@ def via(print):
     warnings.showwarning = oldshowwarning
 
 
+@contextlib.contextmanager
+def ignore(category=Warning):
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", category)
+        yield
+
+
 # vim:sw=4:sts=4:et
diff --git a/tests/test_SI.py b/tests/test_SI.py
index 7c42e22e0..bb030e4dc 100644
--- a/tests/test_SI.py
+++ b/tests/test_SI.py
@@ -145,6 +145,11 @@ def test_sum(self):
         self.assertTrue(numpy.all(numpy.sum(SI.units.kg * numpy.arange(6).reshape(2, 3), 0) == SI.units.kg * numpy.array([3, 5, 7])))
         self.assertTrue(numpy.all(numpy.sum(SI.units.kg * numpy.arange(6).reshape(2, 3), 1) == SI.units.kg * numpy.array([3, 12])))
 
+    def test_prod(self):
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3), axis=0) == numpy.array([4, 10, 18]) * SI.units.m**2))
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3), axis=1) == numpy.array([6, 120]) * SI.units.m**3))
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3)) == 720 * SI.units.m**6))
+
     def test_mean(self):
         self.assertTrue(numpy.all(numpy.mean(SI.units.kg * numpy.arange(6).reshape(2, 3), 0) == SI.units.kg * numpy.array([1.5, 2.5, 3.5])))
         self.assertTrue(numpy.all(numpy.mean(SI.units.kg * numpy.arange(6).reshape(2, 3), 1) == SI.units.kg * numpy.array([1, 4])))
diff --git a/tests/test_basis.py b/tests/test_basis.py
index a6675baf8..121229fc7 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -275,7 +275,7 @@ def setUp(self):
         elif self.variant == 'tensor':
             structured, geom = mesh.rectilinear([numpy.linspace(0, 1, 5-i) for i in range(self.ndims)])
             domain = topology.ConnectedTopology(structured.space, structured.references, structured.transforms, structured.opposites, structured.connectivity)
-            nverts = numpy.product([5-i for i in range(self.ndims)])
+            nverts = numpy.prod([5-i for i in range(self.ndims)])
         elif self.variant == 'simplex':
             numpy.random.seed(0)
             nverts = 20
@@ -299,7 +299,7 @@ def test_ndofs(self):
         elif self.btype == 'discont':
             ndofs_by_ref = {
                 element.getsimplex(1)**self.ndims: (self.degree+1)**self.ndims,
-                element.getsimplex(self.ndims): numpy.product(self.degree+numpy.arange(self.ndims)+1) // numpy.product(numpy.arange(self.ndims)+1)}
+                element.getsimplex(self.ndims): numpy.prod(self.degree+numpy.arange(self.ndims)+1) // numpy.prod(numpy.arange(self.ndims)+1)}
             ndofs = sum(ndofs_by_ref[reference] for reference in self.domain.references)
         elif self.degree == 1:
             ndofs = self.nverts
diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py
index d10e4b389..a3e12649d 100644
--- a/tests/test_evaluable.py
+++ b/tests/test_evaluable.py
@@ -512,8 +512,8 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('arctan-complex', evaluable.arctan, numpy.arctan, ANC(4, 4))
 _check('ln', evaluable.ln, numpy.log, POS(4, 4))
 _check('ln-complex', evaluable.ln, numpy.log, NZC(4, 4))
-_check('product', lambda a: evaluable.product(a, 2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
-_check('product-complex', lambda a: evaluable.product(a, 2), lambda a: numpy.product(a, 2), ANC(4, 3, 4))
+_check('product', lambda a: evaluable.product(a, 2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
+_check('product-complex', lambda a: evaluable.product(a, 2), lambda a: numpy.prod(a, 2), ANC(4, 3, 4))
 _check('sum', lambda a: evaluable.sum(a, 2), lambda a: a.sum(2), ANY(4, 3, 4))
 _check('sum-complex', lambda a: evaluable.sum(a, 2), lambda a: a.sum(2), ANC(4, 3, 4))
 _check('transpose1', lambda a: evaluable.transpose(a, [0, 1, 3, 2]), lambda a: a.transpose([0, 1, 3, 2]), ANY(2, 3, 4, 5))
diff --git a/tests/test_function.py b/tests/test_function.py
index e2f21baa5..ae69d124e 100644
--- a/tests/test_function.py
+++ b/tests/test_function.py
@@ -342,10 +342,10 @@ def _check(name, op, n_op, *args):
 _check('sum-bool', lambda a: numpy.sum(function.Array.cast(a > 0), 2), lambda a: (a > 0).sum(2), ANY(4, 3, 4))
 _check('sum-complex', lambda a: numpy.sum(function.Array.cast(a), 2), lambda a: a.sum(2), ANC(4, 3, 4))
 _check('Array_sum', lambda a: function.Array.cast(a).sum(2), lambda a: a.sum(2), ANY(4, 3, 4))
-_check('product', lambda a: numpy.product(function.Array.cast(a), 2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
-_check('product-bool', lambda a: numpy.product(function.Array.cast(a > 0), 2), lambda a: numpy.product((a > 0), 2), ANY(4, 3, 4))
-_check('product-complex', lambda a: numpy.product(function.Array.cast(a), 2), lambda a: numpy.product(a, 2), ANC(4, 3, 4))
-_check('Array_prod', lambda a: function.Array.cast(a).prod(2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
+_check('product', lambda a: numpy.prod(function.Array.cast(a), 2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
+_check('product-bool', lambda a: numpy.prod(function.Array.cast(a > 0), 2), lambda a: numpy.prod((a > 0), 2), ANY(4, 3, 4))
+_check('product-complex', lambda a: numpy.prod(function.Array.cast(a), 2), lambda a: numpy.prod(a, 2), ANC(4, 3, 4))
+_check('Array_prod', lambda a: function.Array.cast(a).prod(2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
 
 _check('dot', lambda a, b: numpy.dot(a, function.Array.cast(b)), numpy.dot, ANY(1, 2, 5), ANY(3, 5, 4))
 _check('dot-complex', lambda a, b: numpy.dot(a, function.Array.cast(b)), numpy.dot, ANC(1, 2, 5), ANC(3, 5, 4))