diff --git a/hw4-pybind/README.md b/hw4-pybind/README.md
index 9f509922..a191c29e 100644
--- a/hw4-pybind/README.md
+++ b/hw4-pybind/README.md
@@ -1,10 +1,12 @@
 
 # Homework 4
 
 ### Exercise 4:
 Can run the example using the command `./particles 365 1 ../init.csv planet 1` from the `build` folder. From C++ the gravitational constant G must be changed accordingly.
 
 The trajectories can be viewed using Paraview, and we observe trajectories as expected. Except for Mercury.
 
 ### Exercise 5:
 See the file `python_functions.py`.
+Note that we assume 5 integers in the filenames (not 4 as in the trajectory folder).
+We can easily rename multiple files using the command `rename 's/step-0/step-00/g' *` in order to get same naming convention.
diff --git a/hw4-pybind/python_functions.py b/hw4-pybind/python_functions.py
new file mode 100644
index 00000000..6f11d704
--- /dev/null
+++ b/hw4-pybind/python_functions.py
@@ -0,0 +1,34 @@
+import numpy as np
+import os
+
+# readPositions reads the trajectory of a given planet and makes a numpy array out of it
+# returns: 365 × 3 matrix, the columns being the three components of the planet position.
+def readPositions(planet_name, directory):
+    out = np.zeros((365, 3))
+    for s in range(0, 365):
+        filename = directory + "/step-" + str(s).zfill(5) + ".csv"
+        #print(filename)
+        with open(filename) as file:
+            for line in file:
+                line = line.split()
+                for word in line:
+                    #print(word)
+                    if word == planet_name:
+                        out[s, 0] = line[0]
+                        out[s, 1] = line[1]
+                        out[s, 2] = line[2]
+                        break;
+    return out
+
+# computes L2 norm of the error in position
+def computeError(positions, positions_ref):
+    return np.sqrt( np.sum( (positions-positions_ref)**2) )
+
+
+# Testing functions
+planet = "mercury"
+positions     = readPositions(planet, "dumps")
+positions_ref = readPositions(planet, "../trajectories")
+
+error = computeError(positions, positions_ref)
+print(error)