diff --git a/work/week14/particle-pybind/starting_point/generate_input.py b/work/week14/particle-pybind/starting_point/generate_input.py
new file mode 100644
index 0000000..e793655
--- /dev/null
+++ b/work/week14/particle-pybind/starting_point/generate_input.py
@@ -0,0 +1,168 @@
+import os
+import numpy as np
+import math
+import argparse
+import matplotlib.pyplot as plt
+
+
+def isSquare(n):
+    ''' Check if an integer value is a perfect square or not. '''
+    x = n // 2
+    seen = set([x])
+    while x * x != n:
+        x = (x + (n // x)) // 2
+        if x in seen:
+            return False
+        seen.add(x)
+    return True
+
+
+def plotDistributions(hv2d, temp2d):
+    ''' Plot 2D heat and temperature distributions. '''
+    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
+
+    ax = axes[0]
+    ax.set_title('Heat source distribution')
+    ax.set_xlabel('Y')
+    ax.set_ylabel('X')
+    ax.imshow(hv2d, extent=[-1, 1, -1, 1])
+
+    ax = axes[1]
+    ax.set_title('Temperature distribution')
+    ax.set_xlabel('Y')
+    ax.set_ylabel('X')
+    ax.imshow(temp2d, extent=[-1, 1, -1, 1])
+
+    return fig
+
+
+def main():
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-n', '--number', type=int, default=None,
+                        help='Number of particles')
+    parser.add_argument('-o', '--outfile', type=str, default=None,
+                        help='Name of output file')
+    parser.add_argument('-H', '--heatfield', type=str, default='circular',
+                        help='Heat distribution type ("null", "line" or "circular")')
+    parser.add_argument('-T', '--tempfield', type=str, default=None,
+                        help='Initial temperature field ("random" or "homogeneous"),\
+                        only used for non-test heatfields')
+    parser.add_argument('-t', '--test', default=False, action='store_true',
+                        help='Test configuration (temperature field automtcially generated')
+    parser.add_argument('-R', '--radius', type=float, default=1 / 3,
+                        help='Heat source radius (if heatfield is set to "circular")')
+    parser.add_argument('-p', '--plot', default=False, action='store_true',
+                        help='Plot heat and temperature distribution')
+    parser.add_argument('--hr', type=float, default=1,
+                        help='Specific heat rate')
+
+    args = parser.parse_args()
+
+    # Check for input validity
+    test_error_str = 'cannot be specified for test configurations'
+    if args.test:
+        if args.outfile is not None:
+            raise ValueError('output file {}'.format(test_error_str))
+        if args.tempfield is not None:
+            raise ValueError('temperature field {}'.format(test_error_str))
+        if args.number is not None:
+            raise ValueError('number of particles {}'.format(test_error_str))
+
+    # Get number of particles and square root of it
+    N = int(args.number) if args.number is not None else 16**2
+    if not(isSquare(N)):
+        raise ValueError('number of particles must be square')
+    sqrtN = int(math.sqrt(N))
+
+    # Set initial positions on a 2D rectilinear grid (z = 0)
+    x = np.linspace(-1, 1, sqrtN)
+    pos2d = np.array(np.meshgrid(x, x)).T
+
+    # Set heat volumetric source according to user input
+    hv2d = np.zeros((sqrtN, sqrtN))
+    if args.heatfield == 'null':
+        if args.test:
+            # Test configuration: set constant temperature field
+            theta2d = np.ones((sqrtN, sqrtN))
+            outfile = 'testnull.csv'
+
+    elif args.heatfield == 'line':
+        # Set punctual (dirac) heat distribution
+        hv2d[int(sqrtN / 4), :] = -sqrtN
+        hv2d[int(3 * sqrtN / 4), :] = sqrtN
+
+        # Test configuration: set corresponding temperature field from analytical solution
+        if args.test:
+            outfile = 'testline.csv'
+            L = 2
+            freqs = np.fft.fftfreq(sqrtN) * 2 * np.pi / L * sqrtN
+            freqs_2d = np.einsum('i,j->ij', np.ones(sqrtN), freqs)
+            freqs = freqs_2d**2 + freqs_2d.T**2
+            freqs[0, 0] = 1.
+            T = np.fft.fft2(hv2d) / freqs
+            T[0, 0] = 0.
+            T = np.fft.ifft2(T)
+            theta2d = T.real
+
+    elif args.heatfield == 'circular':
+        if args.test:
+            raise ValueError('no test configuration available for "circular" heat field')
+
+        # Set circular heat distribution
+        R = args.radius
+        if R < 0:
+            raise ValueError('radius must be positive')
+        for i in range(sqrtN):
+            for j in range(sqrtN):
+                if (x[i]**2 + x[j]**2 < R):
+                    hv2d[i, j] = sqrtN
+
+    else:
+        raise ValueError('heatfield must be "null", "line" or "circular"')
+
+    if not args.test:
+        # Set temperature field according to user input
+        if args.tempfield == 'random' or args.tempfield is None:
+            # Random field (between zero and 2)
+            theta2d = np.random.random((sqrtN, sqrtN)) * 2
+        elif args.tempfield == 'homogeneous':
+            # Homogeneous field (all equal to one)
+            theta2d = np.ones((sqrtN, sqrtN))
+        else:
+            raise ValueError('tempfield must be "random" or "homogeneous"')
+
+        # Set output file according to user input
+        outfile = args.outfile if args.outfile is not None else 'step-00000.csv'
+
+    # Reshape matrices into serialized vectors
+    positions = np.hstack((pos2d.reshape(N, 2), np.zeros((N, 1))))
+    hv = np.reshape(hv2d, (hv2d.size, 1))
+    theta = np.reshape(theta2d, (theta2d.size, 1))
+
+    # Generate serialized vectors for other, constant quantities
+    hr = args.hr * np.ones((N, 1))
+    velocities = np.zeros((N, 3))
+    forces = np.zeros((N, 3))
+    masses = np.ones((N, 1))
+
+    # Save data into csv file
+    file_data = np.hstack((positions, velocities, forces, masses, theta, hr, hv))
+    outdir = os.path.join(os.getcwd(), 'build', 'dumps')
+    if not os.path.isdir(outdir):
+        os.makedirs(outdir)
+    filepath = os.path.join(os.getcwd(), 'build', 'dumps', outfile)
+    np.savetxt(filepath, file_data, delimiter=" ")
+
+    # Plot the heat and temperature fields
+    if args.plot:
+        plotDistributions(hv2d, theta2d)
+        plt.show()
+
+
+if __name__ == '__main__':
+    try:
+        main()
+    except ValueError as e:
+        print('Error:', e)
+        quit()
diff --git a/work/week14/particle-pybind/starting_point/pypart_pybind.cpp b/work/week14/particle-pybind/starting_point/pypart_pybind.cpp
index e15b3a9..60ffd6f 100644
--- a/work/week14/particle-pybind/starting_point/pypart_pybind.cpp
+++ b/work/week14/particle-pybind/starting_point/pypart_pybind.cpp
@@ -1,84 +1,90 @@
 #include <functional>
 #include <iostream>
 #include <pybind11/pybind11.h>
 
-#include "my_types.hh"
 #include "system_evolution.hh"
-#include "system.hh"
 #include "csv_writer.hh"
 #include "particles_factory_interface.hh"
 #include "ping_pong_balls_factory.hh"
 #include "material_points_factory.hh"
 #include "planets_factory.hh"
 #include "compute_temperature.hh"
 
 namespace py = pybind11;
 
 PYBIND11_MODULE(pypart, m) {
 
-// Wrap the SystemEvolution class
-py::class_<System>(m, "System")
-    .def(py::init<>()) // constructor
-    ;
+    // Bind the SystemEvolution class
+    py::class_<SystemEvolution>(m, "SystemEvolution")
+        // .def(py::init<std::unique_ptr<System>>())
+        .def("addCompute", (&SystemEvolution::addCompute))
+        .def("setNSteps", &SystemEvolution::setNSteps)
+        .def("setDumpFreq", (&SystemEvolution::setDumpFreq))
+        .def("evolve", (&SystemEvolution::setDumpFreq))
+        .def("getSystem", (&SystemEvolution::getSystem))
+        ;
 
-// Wrap the SystemEvolution class
-py::class_<SystemEvolution>(m, "SystemEvolution")
-    // .def(py::init<std::unique_ptr<System>>())
-    .def("addCompute",(&SystemEvolution::addCompute))
-    ;
+    // Bind the ParticlesFactoryInterface class
+    py::class_<ParticlesFactoryInterface>(m, "ParticlesFactoryInterface")
+        // .def(py::init<>()) // constructor
+        .def("getInstance", &ParticlesFactoryInterface::getInstance)
+        .def("createSimulation", (&ParticlesFactoryInterface::createSimulation))
+        ;
 
-// Wrap the ParticlesFactoryInterface class
-py::class_<ParticlesFactoryInterface>(m, "ParticlesFactoryInterface")
-    // .def(py::init<>()) // constructor
-    .def("getInstance", &ParticlesFactoryInterface::getInstance)
-    .def("createSimulation",(&ParticlesFactoryInterface::createSimulation))
-    ;
+    // Bind the PlanetsFactory class
+    py::class_<PlanetsFactory, ParticlesFactoryInterface>(m, "PlanetsFactory")
+        // .def(py::init<>()) // constructor
+        .def("getInstance", &PlanetsFactory::getInstance)
+        .def("createSimulation", &PlanetsFactory::createSimulation)
+        ;
 
-// Wrap the PlanetsFactory class
-py::class_<PlanetsFactory, ParticlesFactoryInterface>(m, "MaterialPointsFactory")
-    // .def(py::init<>()) // constructor
-    .def("getInstance", &PlanetsFactory::getInstance)
-    .def("createSimulation", &PlanetsFactory::createSimulation)
-    ;
+    // Bind the PingPongBallsFactory class
+    py::class_<PingPongBallsFactory, ParticlesFactoryInterface>(m, "PingPongBallsFactory")
+        // .def(py::init<>()) // constructor
+        .def("getInstance", &PingPongBallsFactory::getInstance)
+        .def("createSimulation", &PingPongBallsFactory::createSimulation)
+        ;
 
-// Wrap the PingPongBallsFactory class
-py::class_<PingPongBallsFactory, ParticlesFactoryInterface>(m, "MaterialPointsFactory")
-    // .def(py::init<>()) // constructor
-    .def("getInstance", &PingPongBallsFactory::getInstance)
-    .def("createSimulation", &PingPongBallsFactory::createSimulation)
-    ;
+    // Bind the MaterialPointsFactory class
+    py::class_<MaterialPointsFactory, ParticlesFactoryInterface>(m, "MaterialPointsFactory")
+        // .def(py::init<>()) // constructor
+        .def("getInstance", &MaterialPointsFactory::getInstance)
 
-// Wrap the MaterialPointsFactory class
-py::class_<MaterialPointsFactory, ParticlesFactoryInterface>(m, "MaterialPointsFactory")
-    // .def(py::init<>()) // constructor
-    .def("getInstance", &MaterialPointsFactory::getInstance)
-    .def("createSimulation",
-         (SystemEvolution & (MaterialPointsFactory::*)(const std::string &, Real))
-         &MaterialPointsFactory::createSimulation)  // static casting for "classic" overloaded method
-    .def("createSimulation",
-         py::overload_cast<const std::string &, Real, py::function>
-         (&MaterialPointsFactory::createSimulation<py::function>)) // dynamic casting for "functor" overloaded method
-    ;
+        // static casting for "classic" overloaded method
+        .def("createSimulation",
+             (SystemEvolution & (MaterialPointsFactory::*)(const std::string &, Real))
+             &MaterialPointsFactory::createSimulation)
 
-//Make a shared Compute Method
-py::class_<Compute>(m, "Compute")
-    // .def(py::init<>()) // constructor
-    ;
+        // dynamic casting for "functor" overloaded method
+        .def("createSimulation",
+             py::overload_cast<const std::string &, Real, py::function>
+             (&MaterialPointsFactory::createSimulation<py::function>))
+        ;
 
-//Creates the ComputeTemeprature class and its functions
-py::class_<ComputeTemperature, Compute>(m, "ComputeTemperature")
-    .def(py::init<>()) // constructor)
-    .def_property("conductivity_property", &ComputeTemperature::setConductivity,
-                  &ComputeTemperature::getConductivity) // method only in Dog;
-    .def_property("capacity_property", &ComputeTemperature::setCapacity,
-                  &ComputeTemperature::getCapacity)
-    .def_property("density_property", &ComputeTemperature::setDensity,
-                  &ComputeTemperature::getDensity)
-    .def_property("l_property", &ComputeTemperature::setL,
-                  &ComputeTemperature::getL)
-    .def_property("deltat_property", &ComputeTemperature::setDeltat,
-                  &ComputeTemperature::getDeltat)
-    ;
-}
+    // Bind the generic Compute class (shared compute method ???)
+    py::class_<Compute>(m, "Compute")
+        // .def(py::init<>()) // constructor
+        ;
+
+    // Bind the ComputeTemperature class
+    py::class_<ComputeTemperature, Compute>(m, "ComputeTemperature")
+        .def(py::init<>())  // constructor)
+        .def_property("conductivity_property", &ComputeTemperature::setConductivity,
+                      &ComputeTemperature::getConductivity) // method only in Dog;
+        .def_property("capacity_property", &ComputeTemperature::setCapacity,
+                      &ComputeTemperature::getCapacity)
+        .def_property("density_property", &ComputeTemperature::setDensity,
+                      &ComputeTemperature::getDensity)
+        .def_property("l_property", &ComputeTemperature::setL,
+                      &ComputeTemperature::getL)
+        .def_property("deltat_property", &ComputeTemperature::setDeltat,
+                      &ComputeTemperature::getDeltat)
+        ;
 
+    // Bind the CsvWriter class
+    py::class_<CsvWriter>(m, "CsvWriter")
+        .def(py::init<const std::string&>())  // constructor
+        .def("write", &CsvWriter::write)
+        ;
 
+}