diff --git a/docs/source/graphene_example.md b/docs/source/graphene_example.md
index cd3a2657e5c3263d7af07752dfe0ae01230b1fa2..4d38ec51d921d09a6d66f1a6175dce1224989ff9 100644
--- a/docs/source/graphene_example.md
+++ b/docs/source/graphene_example.md
@@ -66,7 +66,7 @@ After we have defined the guess, we feed it together with the model into the mea
 
 ## Creating a phase diagram of the gap
 
-We can now create a phase diagram of the gap of the interacting solution. In order to calculate the gap we first create a function which takes a hopping dictionary and a Fermi energy and returns the indirect gap. The gap is defined as the difference between the highest occupied and the lowest unoccupied energy level. We will use a dense k-grid to calculate the gap. In order to obtain the Hamiltonian on a dense k-grid, we use the `tb_to_khamvector` function from the `transforms` module.
+We can now create a phase diagram of the gap of the interacting solution. In order to calculate the gap we first create a function which takes a hopping dictionary and a Fermi energy and returns the indirect gap. The gap is defined as the difference between the highest occupied and the lowest unoccupied energy level. We will use a dense k-grid to calculate the gap. In order to obtain the Hamiltonian on a dense k-grid, we use the `tb_to_kgrid` function from pymf.
 
 ```{code-cell} ipython3
 def compute_gap(h, fermi_energy=0, nk=100):
@@ -104,7 +104,7 @@ def compute_phase_diagram(Us, Vs, int_builder, h_0):
   mf_sols = np.asarray(mf_sols).reshape((len(Us), len(Vs)))
   return gaps, mf_sols
 ```
-We chose to initialize a new guess for each $U$ value, but not for each $V$ value. Instead, for consecutive $V$ values we use the previous mean-field solution as a guess. We do this because the mean-field solution is expected to be smooth in the interaction strength and therefore by using an inspired guess we can speed up the calculation.
+We chose to initialize a new guess for each new $U$ and $V$ value. For certain mean-field problems, one might want to reuse the mean-field solution of a nearby parameter as the next guess in order to speed up computations. However, as the size of this system is still small, we can afford to initialize a new guess for each new $U$ and $V$ value.
 
 We can now compute the phase diagram and then plot it
 
diff --git a/docs/source/hubbard_1d.md b/docs/source/hubbard_1d.md
index 9c5b75e54cdf9171ee749cf051da152762f020c4..22b18af9b4c382c073c716b55e34500fed5ac311 100644
--- a/docs/source/hubbard_1d.md
+++ b/docs/source/hubbard_1d.md
@@ -21,12 +21,15 @@ First, the basic imports are done.
 import numpy as np
 import matplotlib.pyplot as plt
 import pymf
+from tqdm import tqdm
 ```
 
 After this, we start by constructing the non-interacting Hamiltonian. As we expect the ground state to be an antiferromagnet, we build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian is then:
+
 $$
 H_0 = \sum_i c_{i, B}^{\dagger}c_{i, A} + c_{i, A}^{\dagger}c_{i+1, B} + h.c.
 $$
+
 We write down the spinful part by simply taking $H_0(k) \otimes \mathbb{1}$.
 
 To translat ethis Hamiltonian into a tight-binding model, we specify the hopping vectors together with the hopping amplitudes. We ensure that the Hamiltonian is hermitian by letting the hopping amplitudes from $A$ to $B$ be the complex conjugate of the hopping amplitudes from $B$ to $A$.
@@ -36,7 +39,7 @@ hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
 h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
 ```
 
-We verify this tight-binding model by plotting the band structure and observing the two bands due to the Brillouin zone folding. In order to do this we transform the tight-binding model into a Hamiltonian on a k-point grid using the `tb_to_khamvector` and then diagonalize it.
+We verify this tight-binding model by plotting the band structure and observing the two bands due to the Brillouin zone folding. In order to do this we transform the tight-binding model into a Hamiltonian on a k-point grid using the `tb_to_kgrid` and then diagonalize it.
 
 ```{code-cell} ipython3
 # Set number of k-points
@@ -70,10 +73,10 @@ This is simply constructed by writing:
 
 ```{code-cell} ipython3
 U=0.5
-h_int = {(0,): U * np.kron(np.ones((2, 2)), np.eye(2)),}
+h_int = {(0,): U * np.kron(np.eye(2), np.ones((2,2))),}
 ```
 
-In order to find a mean-field solution, we combine the non interacting Hamiltonian with the interaction Hamiltonian and the relevant filling into a `Model` object. We then generate a starting guess for the mean-field solution and solve the model using the `solver` function. It is important to note that the guess will influence the possible solutions which the `solver` can find in the mean-field procedure. The `generate_guess` function generates a random Hermitian tight-binding dictionary, with the keys provided as hopping vectors and the values of the size as specified.
+In order to find a mean-field solution, we combine the non interacting Hamiltonian with the interaction Hamiltonian and the relevant filling into a `Model` object. We then generate a starting guess for the mean-field solution and solve the model using the `solver` function. It is important to note that the guess will influence the possible solutions which the `solver` can find in the mean-field procedure. The `generate_guess` function generates a random Hermitian tight-binding dictionary, with the keys provided as hopping vectors and the values of the size as specified. We specifically choose the keys, meaning the hopping vectors, for the `guess` to be the same as the hopping vectors for the interaction Hamiltonian. This is because we do not expect the mean-field solution to contain terms more than the hoppings from the interacting part.
 
 ```{code-cell} ipython3
 filling = 2
@@ -91,12 +94,12 @@ def compute_sol(U, h_0, nk, filling=2):
     }
     guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
     full_model = pymf.Model(h_0, h_int, filling)
-    mf_sol = pymf.solver(full_model, guess, nk=nk, optimizer_kwargs={"M":0})
+    mf_sol = pymf.solver(full_model, guess, nk=nk)
     return pymf.add_tb(h_0, mf_sol)
 
 
 def compute_gap_and_vals(full_sol, nk_dense, fermi_energy=0):
-    h_kgrid = pymf.tb_to_khamvector(full_sol, nk_dense)
+    h_kgrid = pymf.tb_to_kgrid(full_sol, nk_dense)
     vals = np.linalg.eigvalsh(h_kgrid)
 
     emax = np.max(vals[vals <= fermi_energy])
@@ -119,3 +122,20 @@ def compute_phase_diagram(
 
     return np.asarray(gaps, dtype=float), np.asarray(vals)
 ```
+
+Now we run this function over a larger range of $U$ values and plot the gap as a function of $U$.
+
+```{code-cell} ipython3
+Us = np.linspace(0, 8, 40, endpoint=True)
+nk, nk_dense = 40, 100
+gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
+```
+
+```{code-cell} ipython3
+plt.plot(Us, gap, c="k")
+plt.xlabel("$U$")
+plt.ylabel("Gap")
+plt.show()
+```
+
+As expected, there is a critical $U$ around $U=1$ where the gap opens. Furthermore, at large $U$, the gap scales linearly with $U$.