diff --git a/kwant/physics/tests/test_symmetry.py b/kwant/physics/tests/test_symmetry.py
index a2c97f560292f3299b67eed56dd9f986517801d1..79f96f09966aea4bf9eb5532b4c8adb4bd7ec82b 100644
--- a/kwant/physics/tests/test_symmetry.py
+++ b/kwant/physics/tests/test_symmetry.py
@@ -38,7 +38,7 @@ def test_projectors():
     projectors = [sparse.coo_matrix(p) for p in projectors]
     DiscreteSymmetry(projectors=projectors)
     # Break one projector - setting them should now throw an error.
-    projectors = [p.todense() for p in projectors]
+    projectors = [p.toarray() for p in projectors]
     projectors[1][2, :] = 0
     assert not np.allclose(sum(projector.dot(projector.conj().T) for projector
                                in projectors), np.eye(6))
@@ -54,8 +54,8 @@ def test_projectors():
     # Check that set_symmetry removes zero columns properly.
     p0 = p0[:, 1::2]
     p1 = p1[:, ::2]
-    assert np.allclose(symmetry.projectors[0].todense(), p0)
-    assert np.allclose(symmetry.projectors[1].todense(), p1)
+    assert np.allclose(symmetry.projectors[0].toarray(), p0)
+    assert np.allclose(symmetry.projectors[1].toarray(), p1)
 
 def test_set_discrete_symm():
     n = 20
@@ -103,11 +103,11 @@ def test_set_discrete_symm():
 
     # Check that with two symmetries specified, the third one is computed.
     symm = DiscreteSymmetry(time_reversal=t_mat, particle_hole=p_mat)
-    assert np.allclose(symm.chiral.todense(), c_mat.todense())
+    assert np.allclose(symm.chiral.toarray(), c_mat.toarray())
     DiscreteSymmetry(time_reversal=t_mat, chiral=c_mat)
-    assert np.allclose(symm.particle_hole.todense(), p_mat.todense())
+    assert np.allclose(symm.particle_hole.toarray(), p_mat.toarray())
     DiscreteSymmetry(particle_hole=p_mat, chiral=c_mat)
-    assert np.allclose(symm.time_reversal.todense(), t_mat.todense())
+    assert np.allclose(symm.time_reversal.toarray(), t_mat.toarray())
 
 
 def test_projectors_and_symmetry():
diff --git a/kwant/solvers/mumps.py b/kwant/solvers/mumps.py
index 76c53e25e7a07ba091f2a1593e976c1494c5fdfa..7d0dec0be1652ec07bdedb18548a138f4cd1dcba 100644
--- a/kwant/solvers/mumps.py
+++ b/kwant/solvers/mumps.py
@@ -115,7 +115,7 @@ class Solver(common.SparseSolver):
             tmprhs = b[:, j:min(j + self.nrhs, b.shape[1])]
 
             if not self.sparse_rhs:
-                tmprhs = tmprhs.todense()
+                tmprhs = tmprhs.toarray()
             sols.append(solve(tmprhs)[kept_vars, :])
 
         return np.concatenate(sols, axis=1)
diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py
index c43061bdf22503a121f33cc251757330106738fa..b4582f34492edac8bc9322d8dd1532d20ce38537 100644
--- a/kwant/solvers/sparse.py
+++ b/kwant/solvers/sparse.py
@@ -101,7 +101,7 @@ class Solver(common.SparseSolver):
         sols = []
         vec = np.empty(b.shape[0], complex)
         for j in range(b.shape[1]):
-            vec[:] = b[:, j].todense().flatten()
+            vec[:] = b[:, j].toarray().flatten()
             sols.append(factorized_a(vec)[kept_vars])
 
         return np.asarray(sols).transpose()
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 9c83ac594e698fac82f36bdb63e1197acb25035b..14e2b8c2d86bcecd6627fdf9ec59c9199f966d78 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1130,24 +1130,24 @@ def test_discrete_symmetries():
 
     sym = syst.finalized().discrete_symmetry(args=[0])
     for proj, should_be in zip(sym.projectors, np.identity(3)):
-        assert np.allclose(proj.todense(), should_be.reshape((3, 1)))
-    assert np.allclose(sym.time_reversal.todense(), np.identity(3))
+        assert np.allclose(proj.toarray(), should_be.reshape((3, 1)))
+    assert np.allclose(sym.time_reversal.toarray(), np.identity(3))
     syst.conservation_law = lambda site, p: cons_law[site.family]
     sym = syst.finalized().discrete_symmetry(args=[0])
     for proj, should_be in zip(sym.projectors, np.identity(3)):
-        assert np.allclose(proj.todense(), should_be.reshape((-1, 1)))
+        assert np.allclose(proj.toarray(), should_be.reshape((-1, 1)))
 
     syst = builder.Builder(conservation_law=np.diag([-1, 1]))
     syst[lat(1)] = np.identity(2)
     sym = syst.finalized().discrete_symmetry()
     for proj, should_be in zip(sym.projectors, np.identity(2)):
-        assert np.allclose(proj.todense(), should_be.reshape((-1, 1)))
+        assert np.allclose(proj.toarray(), should_be.reshape((-1, 1)))
 
     syst = builder.Builder(conservation_law=1)
     syst[lat2(1)] = 0
     sym = syst.finalized().discrete_symmetry()
     [proj] = sym.projectors
-    assert np.allclose(proj.todense(), [[1]])
+    assert np.allclose(proj.toarray(), [[1]])
 
     syst = kwant.Builder(conservation_law=np.diag([-1, 1, -1, 1]))
 
@@ -1156,17 +1156,17 @@ def test_discrete_symmetries():
     sym = syst.finalized().discrete_symmetry()
     p1 = np.zeros((4, 2))
     p1[0, 0] = p1[2, 1] = 1
-    assert np.allclose(sym.projectors[0].todense(), p1)
+    assert np.allclose(sym.projectors[0].toarray(), p1)
     p2 = np.zeros((4, 2))
     p2[1, 0] = p2[3, 1] = 1
-    assert np.allclose(sym.projectors[1].todense(), p2)
+    assert np.allclose(sym.projectors[1].toarray(), p2)
 
     # test parameter passing to conservation_law
     syst = builder.Builder(conservation_law=lambda site, b: b)
     syst[lat2(1)] = 0
     sym = syst.finalized().discrete_symmetry(params=dict(a=None, b=1))
     [proj] = sym.projectors
-    assert np.allclose(proj.todense(), [[1]])
+    assert np.allclose(proj.toarray(), [[1]])
 
 
 def test_argument_passing():
diff --git a/kwant/tests/test_operator.py b/kwant/tests/test_operator.py
index 659d5d772f82f2353c5ef8b205dd4375cd846501..e38e1177c7620fab0ccd4ca69747136bd635b74b 100644
--- a/kwant/tests/test_operator.py
+++ b/kwant/tests/test_operator.py
@@ -387,13 +387,13 @@ def test_tocoo():
     assert isinstance(op.tocoo(), coo_matrix)
 
     # Constant and non-constant values.
-    assert np.all(op.tocoo().todense() == np.eye(2))
+    assert np.all(op.tocoo().toarray() == np.eye(2))
     op = ops.Density(syst, lambda site: 1)
-    assert np.all(op.tocoo().todense() == np.eye(2))
+    assert np.all(op.tocoo().toarray() == np.eye(2))
 
     # Correct treatment of where
     op = ops.Density(syst, where=[lat1(0)])
-    assert np.all(op.tocoo().todense() == [[1, 0], [0, 0]])
+    assert np.all(op.tocoo().toarray() == [[1, 0], [0, 0]])
 
     # No accidental transpose.
     syst = kwant.Builder()
@@ -401,7 +401,7 @@ def test_tocoo():
     syst[lat2(0)] = lambda site, paramerer: np.eye(2)
     syst = syst.finalized()
     op = ops.Density(syst, [[1, 1], [0, 1]], check_hermiticity=False)
-    assert np.all(op.tocoo().todense() == [[1, 1], [0, 1]])
+    assert np.all(op.tocoo().toarray() == [[1, 1], [0, 1]])
 
     op = ops.Density(syst, lambda site, p: [[1, 1], [0, 1]],
                      check_hermiticity=False)
@@ -427,7 +427,7 @@ def test_arg_passing(A):
     act_should_be = op.act(wf, args=canonical_args)
     has_tocoo = hasattr(op, 'tocoo')
     if has_tocoo:
-        tocoo_should_be = op.tocoo(args=canonical_args).todense()
+        tocoo_should_be = op.tocoo(args=canonical_args).toarray()
 
     with raises(TypeError) as exc:
         op(wf, args=canonical_args, params=params)
@@ -449,7 +449,7 @@ def test_arg_passing(A):
         act_should_be, op.act(wf, params=params))
     if has_tocoo:
         np.testing.assert_array_equal(
-            tocoo_should_be, op.tocoo(params=params).todense())
+            tocoo_should_be, op.tocoo(params=params).toarray())
     # after binding
     op2 = op.bind(params=params)
     np.testing.assert_array_equal(
@@ -458,7 +458,7 @@ def test_arg_passing(A):
         act_should_be, op2.act(wf))
     if has_tocoo:
         np.testing.assert_array_equal(
-            tocoo_should_be, op2.tocoo().todense())
+            tocoo_should_be, op2.tocoo().toarray())
 
     # system and onsite having different args
     def onsite(site, flip):
@@ -477,7 +477,7 @@ def test_arg_passing(A):
         act_should_be, op.act(wf, params=params))
     if has_tocoo:
         np.testing.assert_array_equal(
-            tocoo_should_be, op.tocoo(params=params).todense())
+            tocoo_should_be, op.tocoo(params=params).toarray())
     # after binding
     op2 = op.bind(params=params)
     np.testing.assert_array_equal(
@@ -486,7 +486,7 @@ def test_arg_passing(A):
         act_should_be, op2.act(wf))
     if has_tocoo:
         np.testing.assert_array_equal(
-            tocoo_should_be, op2.tocoo().todense())
+            tocoo_should_be, op2.tocoo().toarray())
 
 
 def random_onsite(i):
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 76ac072d6c7ec9d6e18b7fbafaea01a7eca588e4..9f5216383b65ed778abad5841fe1f0a3218ae8e5 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -36,7 +36,7 @@ def test_hamiltonian_submatrix():
 
     mat = syst2.hamiltonian_submatrix(sparse=True)
     assert sparse.isspmatrix_coo(mat)
-    mat = mat.todense()
+    mat = mat.toarray()
     mat = mat[perm, :]
     mat = mat[:, perm]
     np.testing.assert_array_equal(mat, mat_should_be)
@@ -45,7 +45,7 @@ def test_hamiltonian_submatrix():
     np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
 
     mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]], sparse=True)
-    mat = mat.todense()
+    mat = mat.toarray()
     np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
 
     # Test for correct treatment of matrix input.
@@ -57,7 +57,7 @@ def test_hamiltonian_submatrix():
     syst[chain(2), chain(1)] = np.array([[3j]])
     syst2 = syst.finalized()
     mat_dense = syst2.hamiltonian_submatrix()
-    mat_sp = syst2.hamiltonian_submatrix(sparse=True).todense()
+    mat_sp = syst2.hamiltonian_submatrix(sparse=True).toarray()
     np.testing.assert_array_equal(mat_sp, mat_dense)
 
     # Test precalculation of modes.