diff --git a/docs/examples/ex41.py b/docs/examples/ex41.py
index d8fca22a..385a4cf8 100644
--- a/docs/examples/ex41.py
+++ b/docs/examples/ex41.py
@@ -19,6 +19,7 @@
 on mixed meshes is work-in-progress.
 
 """
+import numpy as np
 from pathlib import Path
 from skfem import *
 from skfem.models import laplace, unit_load
@@ -35,7 +36,7 @@
 A = asm(laplace, basis)
 f = asm(unit_load, basis)
 
-y = solve(*condense(A, f, D=out[0]['boundary']['line']).astype(np.int64))
+y = solve(*condense(A, f, D=out[0]['boundary']['line'].astype(np.int64)))
 
 def visualize():
     from skfem.visuals.matplotlib import plot, draw
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index d67d099c..4b1d6e0f 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -588,7 +588,16 @@ def _remove_duplicate_nodes(p, t):
         tmp = np.ascontiguousarray(p.T)
         tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1]),
                                   return_index=True, return_inverse=True)
-        return p[:, ixa], ixb[t]
+        return p[:, ixa], Mesh._squeeze_if(ixb[t])
+
+    @staticmethod
+    def _squeeze_if(arr):
+        # Workaround for the additional dimension introduced in
+        # numpy 2.0 for the output of np.unique when using
+        # return_index=True
+        if len(arr.shape) > 2:
+            return arr.squeeze(axis=2)
+        return arr
 
     def __iter__(self):
         return iter((self.doflocs, self.t))
@@ -609,8 +618,8 @@ def __matmul__(self, other):
                                     return_index=True, return_inverse=True)
             p = p[:, ixa]
             return [
-                cls(p, ixb[self.t]),
-                *[type(m)(p, ixb[m.t + self.p.shape[1]])
+                cls(p, self._squeeze_if(ixb[self.t])),
+                *[type(m)(p, self._squeeze_if(ixb[m.t + self.p.shape[1]]))
                   for i, m in enumerate(other)],
             ]
         raise NotImplementedError