diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index dc3e665dd3ef9505cacb852226943ca946e82e62..396b9175e8a59a4e5f14169036710d026d3f1048 100644
--- a/kwant/_system.pyx
+++ b/kwant/_system.pyx
@@ -11,6 +11,7 @@ import tinyarray as ta
 import numpy as np
 from scipy import sparse as sp
 from itertools import chain
+import types
 
 from .graph.core cimport CGraph, gintArraySlice
 from .graph.defs cimport gint
@@ -327,3 +328,13 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
         mat = func(ham, args, self.graph, diag, from_sites,
                    n_by_to_site, to_norb, to_off, from_norb, from_off)
     return (mat, to_norb, from_norb) if return_norb else mat
+
+# workaround for Cython functions not having __get__ and
+# Python 3 getting rid of unbound methods
+cdef class HamiltonianSubmatrix:
+
+    def __get__(self, obj, objtype):
+        if obj is None:
+            return hamiltonian_submatrix
+        else:
+            return types.MethodType(hamiltonian_submatrix, obj)
diff --git a/kwant/system.py b/kwant/system.py
index ec68b9bf797c34554458e8244fa7ae90f2e2aac8..45d0b473395b2a189ecf60de62078e815ebca4e8 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -48,8 +48,7 @@ class System(object, metaclass=abc.ABCMeta):
         pass
 
 # Add a C-implemented function as an unbound method to class System.
-System.hamiltonian_submatrix = types.MethodType(
-    _system.hamiltonian_submatrix, None, System)
+System.hamiltonian_submatrix = _system.HamiltonianSubmatrix()
 
 
 class FiniteSystem(System, metaclass=abc.ABCMeta):