diff --git a/codes/test.py b/codes/kwant_examples.py
similarity index 58%
rename from codes/test.py
rename to codes/kwant_examples.py
index 7e6041637d9a21922f67b6793cdb41b841ba3673..0b25865e55f55992fc885df5751acacad6c84ea4 100644
--- a/codes/test.py
+++ b/codes/kwant_examples.py
@@ -1,47 +1,39 @@
 import kwant
 import numpy as np
-
+from . import utils
 
 s0 = np.identity(2)
 sz = np.diag([1, -1])
 
-def graphene_onsite(U, N_ks):
+def graphene_extended_hubbard():
+    # Create graphene lattice
     graphene = kwant.lattice.general(
-        [[1, 0], [1 / 2, np.sqrt(3) / 2]], [[0, 0], [0, 1 / np.sqrt(3)]]
+        [[1, 0], [1 / 2, np.sqrt(3) / 2]], [[0, 0], [0, 1 / np.sqrt(3)]],
+        norbs=2
     )
     a, b = graphene.sublattices
 
-    # create bulk system
+    # Create bulk system
     bulk_graphene = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))
-    # add sublattice potential
-    m0 = 0
-    bulk_graphene[a.shape((lambda pos: True), (0, 0))] = m0 * sz
-    bulk_graphene[b.shape((lambda pos: True), (0, 0))] = -m0 * sz
-    # add hoppings between sublattices
+    # Set onsite energy to zero
+    bulk_graphene[a.shape((lambda pos: True), (0, 0))] = 0 * s0
+    bulk_graphene[b.shape((lambda pos: True), (0, 0))] = 0 * s0
+    # Add hoppings between sublattices
     bulk_graphene[graphene.neighbors(1)] = s0
 
-    # use kwant wraparound to sample bulk k-space
-    wrapped_syst = kwant.wraparound.wraparound(bulk_graphene).finalized()
-
-    # return a hamiltonian for a given kx, ky
-    @np.vectorize
-    def hamiltonian_return(kx, ky, params={}):
-        ham = wrapped_syst.hamiltonian_submatrix(params={**params, **dict(k_x=kx, k_y=ky)})
-        return ham
+    def onsite_int(site, U):
+        return U * np.ones((2, 2))
+    def nn_int(site1, site2, V):
+        return V * np.ones((2, 2))
 
-    N_k_axis = np.linspace(0, 2 * np.pi, N_ks, endpoint=False)
-    hamiltonians_0 = np.array(
-        [[hamiltonian_return(kx, ky) for kx in N_k_axis] for ky in N_k_axis]
-    )
-
-    # we need block diagonal structure here since opposite spins interact on the same sublattice
-    v_int = U * np.block(
-        [[np.ones((2, 2)), np.zeros((2, 2))], [np.zeros((2, 2)), np.ones((2, 2))]]
+    syst_V = utils.build_interacting_syst(
+        syst=bulk_graphene,
+        lattice = graphene,
+        func_onsite = onsite_int,
+        func_hop = nn_int,
     )
-    # repeat the matrix on a k-grid
-    V = np.array([[v_int for i in range(N_ks)] for j in range(N_ks)])
     
-    return hamiltonians_0, V
+    return bulk_graphene, syst_V
 
 def hubbard_2D(U, N_ks):
     square = kwant.lattice.square(a=1, norbs=2)
@@ -75,12 +67,10 @@ def hubbard_1D(U, N_ks):
     chain = kwant.lattice.chain(a=1, norbs=2)
     # create bulk system
     bulk_hubbard = kwant.Builder(kwant.TranslationalSymmetry(*chain.prim_vecs))
-    bulk_hubbard[chain.shape((lambda pos: True), (0,))] = 0 * np.eye(2)
+    bulk_hubbard[chain.shape((lambda pos: True), (0,))] = 0 * s0
     # add hoppings between lattice points
     bulk_hubbard[chain.neighbors()] = -1
 
-    # use kwant wraparound to sample bulk k-space
-    wrapped_fsyst = kwant.wraparound.wraparound(bulk_hubbard).finalized()
 
     # return a hamiltonian for a given kx, ky
     @np.vectorize
@@ -93,7 +83,4 @@ def hubbard_1D(U, N_ks):
         [hamiltonian_return(kx) for kx in N_k_axis]
         )
 
-    # onsite interactions
-    v_int = U * np.ones((2,2))
-    V = np.array([v_int for i in range(N_ks)])
     return hamiltonians_0, V