diff --git a/Homework4/src/error_minimization.py b/Homework4/src/error_minimization.py
index f011d2b..9b2f7e5 100644
--- a/Homework4/src/error_minimization.py
+++ b/Homework4/src/error_minimization.py
@@ -1,105 +1,135 @@
 import numpy as np
 import argparse
 import matplotlib.pyplot as plt
 import glob
 from scipy import optimize
 
 import main
 
-# Functions
+# Python function
+#################
 
+# Reads the position of a given planet over time and stores it in an array
 def readPosition(planet_name, directory):
 
     file_list = sorted(glob.glob(directory + '/*.csv'))
     m = len(file_list)
     n = 3
     positions = np.zeros([m,n])
     i = 0
     for file_path in file_list:
         data = np.genfromtxt(file_path, delimiter = ' ', dtype=np.str)
         loc = np.where(data == planet_name)
         line = loc[0]
         positions[i, 0] = float(data[line, 0])
         positions[i, 1] = float(data[line, 1])
         positions[i, 2] = float(data[line, 2])
         i = i + 1
 
     return positions
 
+# Computes the l2-norm error between the reference and computed positions of a given planet
 def computeError(positions, positions_ref):
 
     e = 0.0
-    print(len(positions))
     for i in range(len(positions)):
         diff = (positions_ref[i, 0] - positions[i, 0])**2 + (positions_ref[i, 1] - positions[i, 1])**2 + (positions_ref[i, 2] - positions[i, 2])**2
         e = e + diff
     error = np.sqrt(e)
+    
     return error
 
+# Modifies the input file by scaling the velocity components of a given planet
 def generateInput(scale, planet_name, input_filename, output_filename):
 
     data = np.genfromtxt(input_filename, delimiter = ' ', dtype=np.str)
     loc = np.where(data == planet_name)
     line = loc[0]
     data[line, 3] = float(data[line, 3])*scale
     data[line, 4] = float(data[line, 4])*scale
     data[line, 5] = float(data[line, 5])*scale
     np.savetxt(output_filename, data, delimiter=' ', fmt="%s")
 
     return 0
 
+# Calls main.py to run the particle code based on a specified input file
 def launchParticles(input_file, nb_steps, freq):
+    
     particle_type = 'planet'
     timestep = 1
     main.main(nb_steps, freq, input_file, particle_type, timestep)
 
     return 0
 
+# Runs the particle code and computes the l2-norm error based on a specified input file
 def runAndComputeError(scale, *args):
+    
     planet_name = args[0]
     input_file = args[1]
     nb_steps = args[2]
     freq = args[3]
     generateInput(scale, planet_name, input_file, 'scaled_input.csv')
     launchParticles('scaled_input.csv', nb_steps, freq)
     positions = readPosition(planet_name, dir_comp)
     positions_ref = readPosition(planet_name, dir_ref)
     error = computeError(positions, positions_ref)
-    print(error)
 
     return error
 
 # Arguments parser
+##################
+    
 parser = argparse.ArgumentParser(description='Run particles code and computes simulation error')
 parser.add_argument('directory_comp', type=str, help='path to the directory of the computed data')
 parser.add_argument('directory_ref', type=str, help='path to the directory of the reference data')
 parser.add_argument('planet', type=str, help='Name of the planet to be investigated')
 parser.add_argument('input_file', type=str, help='Name of the inputs file')
 parser.add_argument('scale', type=float, help='Initial velocity scaling factor')
 parser.add_argument('nb_steps', type=int, help='Number of time steps')
 parser.add_argument('freq', type=int, help='Outputs dumping frequency')
 
 args = parser.parse_args()
 
 dir_comp = args.directory_comp
 dir_ref = args.directory_ref
 planet_name = args.planet
 input_file = args.input_file
 scale = args.scale
 nb_steps = args.nb_steps
 freq = args.freq
 
 # Minimization
-#positions = readPosition(planet_name, dir_comp)
-#positions_ref = readPosition(planet_name, dir_ref)
-#errors = computeError(positions, positions_ref)
-#print(errors)
+##############
+
+# Stores scaling factor and error value after each minimization iteration
 def reporter(scale_r):
+    
     global scale_int
+    global error_int
     scale_int = np.vstack([scale_int, scale_r])
+    error_r = runAndComputeError(scale_r, planet_name, input_file, nb_steps, freq)
+    error_int = np.vstack([error_int, error_r])
+    
     return 0
 
 scale_0 = scale
 scale_int = scale_0
-solution = optimize.fmin(runAndComputeError, scale_0, args = (planet_name, input_file, nb_steps, freq), callback = reporter, xtol = 0.0001, ftol = 0.0001, maxiter = 10)
-print(scale_int)
\ No newline at end of file
+error_int = runAndComputeError(scale_0, planet_name, input_file, nb_steps, freq)
+solution = optimize.fmin(runAndComputeError, scale_0, args = (planet_name, input_file, nb_steps, freq), callback = reporter, xtol = 0.0001, ftol = 0.0001, maxiter = 100, full_output = True)
+
+# Post-Processing
+#################
+
+fig = plt.figure(1)
+ax = fig.add_subplot(111)
+ax.plot(scale_int, error_int, 'r--o')
+
+# Axis labels
+ax.set_xlabel('scaling', fontsize = 14)
+ax.set_ylabel('error: L-2 norm', fontsize = 14)
+
+# Title of the window
+ax.set_title('%s trajectory error minimization' % planet_name, fontsize = 14)
+
+plt.savefig('minimization_evolution.png')
+plt.show()
\ No newline at end of file