diff --git a/hw3-heat-fft/test_temperature.cc b/hw3-heat-fft/test_temperature.cc
new file mode 100644
index 00000000..44898382
--- /dev/null
+++ b/hw3-heat-fft/test_temperature.cc
@@ -0,0 +1,110 @@
+#include "compute_temperature.hh"
+#include "csv_reader.hh"
+#include "csv_writer.hh"
+#include "material_point.hh"
+#include "material_points_factory.hh"
+#include "system.hh"
+#include <gtest/gtest.h>
+#include <vector>
+#include "vector.hh"
+
+
+/*****************************************************************/
+// Fixture class
+class HomogeneousTemperatures : public ::testing::Test {
+protected:
+    void SetUp() override {
+        MaterialPointsFactory::getInstance();
+        std::vector <MaterialPoint> materialpoints;
+        n_materialpoints = 10;
+        for (UInt i = 0; i < n_materialpoints; ++i) {
+            MaterialPoint p;
+
+            Vector pos[2];
+            pos[0] = n_materialpoints / i;
+            pos[1] = i % n_materialpoints;
+            p.getPosition() = pos[2];
+            p.getTemperature() = 1.;
+            p.getHeatRate() = 0.;
+            materialpoints.push_back(p);
+
+            std::cout << p << std::endl;
+            system.addParticle(std::make_shared<MaterialPoint>(p));
+        };
+    };
+
+    System system;
+    UInt n_materialpoints;
+};
+
+/*****************************************************************/
+// Fixture class
+class VolumetricHeatSource : public ::testing::Test {
+protected:
+    void SetUp() override {
+        MaterialPointsFactory::getInstance();
+        std::vector <MaterialPoint> materialpoints;
+        n_materialpoints = 10;
+        for (UInt i = 0; i < n_materialpoints; ++i) {
+            MaterialPoint p;
+            Vector pos[2];
+
+            x = (n_materialpoints / i) * dx - L/2;
+            //pos[0] = (n_materialpoints / i) * dx - L/2;
+            //pos[1] = (i % n_materialpoints) * dx - L/2;
+            //p.getPosition() = pos[2];
+
+            if (x > 0.5 - 0.005 and x < 0.5 + 0.005) {
+                p.getTemperature() = x;
+                p.getHeatRate() = 1.;
+            } else if (x > -0.5 - 0.005 and x < -0.5 + 0.005) {
+                p.getTemperature() = x;
+                p.getHeatRate() = -1.;
+            } else if (x <= -0.5 - 0.005) {
+                p.getTemperature() = -x-1;
+                p.getHeatRate() = 0.;
+            } else if (x >= +0.5 + 0.005){
+                p.getTemperature() = -x+1;
+                p.getHeatRate() = 0.;
+            } else {
+                p.getTemperature() = x;
+                p.getHeatRate() = 0.;
+            };
+
+            materialpoints.push_back(p);
+
+            std::cout << p << std::endl;
+            system.addParticle(std::make_shared<MaterialPoint>(p));
+        };
+    };
+
+    System system;
+    UInt n_materialpoints;
+    UInt L = 2.0;
+    UInt dx = (std::sqrt(n_materialpoints)-1)/L;
+    UInt x = 0.;
+};
+
+
+/*****************************************************************/
+TEST_F(HomogeneousTemperatures, csv) {
+  CsvWriter writer("tmp_file");
+  writer.compute(system);
+
+  //CsvReader reader("tmp_file");
+  //System s_fromfile;
+  //reader.compute(s_fromfile);
+  //verlet->setDeltaT(5e-3);
+
+  //UInt read_n_materialpoints = s_fromfile.getNbParticles();
+  //EXPECT_EQ(n_materialpoints, s_fromfile.getNbParticles());
+
+  for (UInt i = 0; i < n_materialpoints; ++i) {
+    auto& part1 = system.getParticle(i);
+    auto& pos1 = part1.getPosition();
+    for (UInt j = 0; j < 3; ++j) {
+      ASSERT_NEAR(pos1[j], 0, 1e-12);
+    }
+  }
+}
+/*****************************************************************/