diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 599d46e..92bbed6 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/docs/source/chapters/chapter1.rst
@@ -76,6 +76,8 @@ Copy the following lines into *potentials.py*:
 
 .. code-block:: python
 
+    import numpy as np
+
     def potentials(potential_type, epsilon, sigma, r, derivative=False):
         if potential_type == "Lennard-Jones":
             if derivative:
@@ -84,9 +86,11 @@ Copy the following lines into *potentials.py*:
                 return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
         elif potential_type == "Hard-Sphere":
             if derivative:
-                raise ValueError("Derivative is not defined for Hard-Sphere potential.")
+                # Derivative is not defined for Hard-Sphere potential.
+                # --> return 0
+                return np.zeros(len(r))
             else:
-                return 1000 if r < sigma else 0
+                return np.where(r > sigma, 0, 1000)
         else:
             raise ValueError(f"Unknown potential type: {potential_type}")
 
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 3599d03..e460b16 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -146,9 +146,9 @@ following *run()* method to the *MinimizeEnergy* class:
             # First, meevaluate the initial energy and max force
             self.update_neighbor_lists() # Rebuild neighbor list, if necessary
             self.update_cross_coefficients() # Recalculate the cross coefficients, if necessary
-            if self.step == 0: # At the first step, Epot/MaxF do not exists yet, calculate them both
-                init_Epot = self.compute_potential()
-                forces, init_MaxF = self.compute_force()
+            # Compute Epot/MaxF/force
+            init_Epot = self.compute_potential()
+            forces, init_MaxF = self.compute_force()
             # Save the current atom positions
             init_positions = copy.deepcopy(self.atoms_positions)
             # Move the atoms in the opposite direction of the maximum force
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index e285f7c..b15b715 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -143,7 +143,7 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
                 f.write("ITEM: NUMBER OF ATOMS\n")
                 f.write(str(code.total_number_atoms) + "\n")
                 f.write("ITEM: BOX BOUNDS pp pp pp\n")
-                for dim in np.arange(code.dimensions):
+                for dim in np.arange(3):
                     f.write(str(box_boundaries[dim][0]) + " "
                             + str(box_boundaries[dim][1]) + "\n")
                 cpt = 1