diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 53f2379cd53f17588dc9e59cb0300a11e90c54c3..f34642bb99288d4f53f89098464ca589193f458c 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -15,6 +15,7 @@ import numpy as np
 import tinyarray as ta
 
 import sympy
+from sympy.matrices.matrices import MatrixBase
 from sympy.utilities.lambdify import lambdastr
 from sympy.printing.lambdarepr import LambdaPrinter
 from sympy.printing.precedence import precedence
@@ -211,7 +212,7 @@ def discretize_symbolic(hamiltonian, coords=None, *, locals=None):
 
     onsite_zeros = (0,) * len(coords)
 
-    if not isinstance(hamiltonian, sympy.matrices.MatrixBase):
+    if not isinstance(hamiltonian, MatrixBase):
         hamiltonian = sympy.Matrix([hamiltonian])
         _input_format = 'expression'
     else:
@@ -574,7 +575,7 @@ def _return_string(expr, coords):
 
     expr = expr.subs(map_func_calls)
 
-    if isinstance(expr, sympy.matrices.MatrixBase):
+    if isinstance(expr, MatrixBase):
         # express matrix return values in terms of sums of known matrices,
         # which will be assigned to '_cache_n' in the function body.
         mons = monomials(expr, expr.atoms(sympy.Symbol))
diff --git a/kwant/plotter.py b/kwant/plotter.py
index d3bb4493bc32dca040cdf5352a653adc6a3bea4a..e8e510d3ccc733fdaceb8cda33ccc395461e6a79 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -1498,8 +1498,8 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
             h_p = np.atleast_2d(bound_ham(**p))
             spectrum.append(np.linalg.eigvalsh(h_p))
     # massage masked grid points into a list of NaNs of the appropriate length
-    n_eigvals = len(next(filter(lambda s: s is not None, spectrum)))
-    nan_list = [np.nan] * n_eigvals
+    shape_eigvals = next(filter(lambda s: s is not None, spectrum)).shape
+    nan_list = np.full(shape_eigvals, np.nan)
     spectrum = [nan_list if s is None else s for s in spectrum]
     # make into a numpy array and reshape
     new_shape = [len(v) for v in array_values] + [-1]
@@ -1542,9 +1542,11 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
         # plot_surface cannot directly handle rank-3 values, so we
         # explicitly loop over the last axis
         grid = np.meshgrid(*array_values)
-        for i in range(spectrum.shape[-1]):
-            spec = spectrum[:, :, i].transpose()  # row-major to x-y ordering
-            ax.plot_surface(*(grid + [spec]), cstride=1, rstride=1)
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore', message='Z contains NaN values')
+            for i in range(spectrum.shape[-1]):
+                spec = spectrum[:, :, i].transpose()  # row-major to x-y ordering
+                ax.plot_surface(*(grid + [spec]), cstride=1, rstride=1)
 
     _maybe_output_fig(fig, file=file, show=show)
 
@@ -1687,8 +1689,13 @@ def _interpolate_field(dim, elements, discrete_field, bbox, width,
 
         # Coordinates of the grid points that are within range of the current
         # hopping.
-        coords = np.meshgrid(*[region[d][field_slice[d]] for d in range(dim)],
-                             sparse=True, indexing='ij')
+        coords = np.array(
+            np.meshgrid(
+                *[region[d][field_slice[d]] for d in range(dim)],
+                sparse=True, indexing='ij'
+            ),
+            dtype=object
+        )
 
         # Convert "coords" into scaled distances from pos_offset
         coords -= pos_offsets[i]
diff --git a/kwant/qsymm.py b/kwant/qsymm.py
index a6dae1c4ff1c3183f829ab0fe2dd643c3fbab1d8..1bd12a2ec04c5f6a7c65bfe4feac705efebd0c3c 100644
--- a/kwant/qsymm.py
+++ b/kwant/qsymm.py
@@ -18,6 +18,7 @@ import scipy.linalg as la
 
 try:
     import sympy
+    import sympy.matrices.matrices
     import qsymm
     from qsymm.model import Model, BlochModel, BlochCoeff
     from qsymm.groups import PointGroupElement, ContinuousGroupGenerator
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index ecf33f80003c6aee65ea5c90c3c75b0f47bc469b..9bcfeed072213d96162bc6e0cd9b18ac7e6eed2a 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -814,7 +814,7 @@ class SMatrix(BlockResult):
                 block_offsets.append(block_offset)
         # Symmetry block offsets for all leads - or None if lead does not have
         # blocks.
-        self.block_offsets = block_offsets
+        self.block_offsets = np.array(block_offsets, dtype=object)
         # Pick out symmetry block offsets for in and out leads
         self.in_block_offsets = \
                 np.array(self.block_offsets)[list(self.in_leads)]