diff --git a/docs/source/graphene_example.md b/docs/source/graphene_example.md
index 1931a9ae2322b28c3eb4c2329dc5f752f3eee15c..cd3a2657e5c3263d7af07752dfe0ae01230b1fa2 100644
--- a/docs/source/graphene_example.md
+++ b/docs/source/graphene_example.md
@@ -70,7 +70,7 @@ We can now create a phase diagram of the gap of the interacting solution. In ord
 
 ```{code-cell} ipython3
 def compute_gap(h, fermi_energy=0, nk=100):
-    kham = pymf.tb_to_khamvector(h, nk, ks=None)
+    kham = pymf.tb_to_kgrid(h, nk)
     vals = np.linalg.eigvalsh(kham)
 
     emax = np.max(vals[vals <= fermi_energy])
diff --git a/docs/source/hubbard_1d.md b/docs/source/hubbard_1d.md
index 8f1a7ea0f8989e4f18aa7d79ff83e9281b9b2469..9c5b75e54cdf9171ee749cf051da152762f020c4 100644
--- a/docs/source/hubbard_1d.md
+++ b/docs/source/hubbard_1d.md
@@ -42,7 +42,7 @@ We verify this tight-binding model by plotting the band structure and observing
 # Set number of k-points
 nk = 100
 ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
-hamiltonians_0 = pymf.tb_to_khamvector(h_0, nk, ks=ks)
+hamiltonians_0 = pymf.tb_to_kgrid(h_0, nk)
 
 vals, vecs = np.linalg.eigh(hamiltonians_0)
 plt.plot(ks, vals, c="k")