diff --git a/exercice_3/src/compute_temperature.cc b/exercice_3/src/compute_temperature.cc
index 1ae54ed..8a4b26e 100644
--- a/exercice_3/src/compute_temperature.cc
+++ b/exercice_3/src/compute_temperature.cc
@@ -1,63 +1,63 @@
 #include "compute_temperature.hh"
 #include "fft.hh"
 #include "material_point.hh"
 #include <cmath>
 
 /* -------------------------------------------------------------------------- */
 
 void ComputeTemperature::compute(System &system) {
   // create temperature field matrix --> skech of operations
 
   // casting the matrices
   // preparing the temp and heatrate vectors
   std::vector<Real> temperatures;
   std::vector<Real> heat_rates;
   // getting the vectors
   for (auto &npar : system) {
     auto &part = static_cast<MaterialPoint &>(npar);
     temperatures.push_back(part.getTemperature());
     heat_rates.push_back(part.getHeatRate());
   }
   // creating the matrices out of the vectors
   UInt nb_particles = system.getNbParticles();
   UInt size = std::sqrt(nb_particles);
 
   auto temperatures_matrix = makeMatrix(temperatures);
   auto hv_matrix = makeMatrix(heat_rates);
 
   // step 1: forward fft
   auto temperaturesfft_matrix = FFT::transform(temperatures_matrix);
   auto hvfft_matrix = FFT::transform(hv_matrix);
-  auto qx2qy2 = FFT::computeFrequenciesqx2qy2(size);
 
   // step 2: dteta^/dt = 1/rho*C * (hv_fft - k*Te_fft(qx**2 + qy**2))
+  auto qx2qy2 = FFT::computeFrequenciesqx2qy2(size);
   auto kTe1 = matrix_mult_sm(k, temperaturesfft_matrix);
   auto kTe1qx2qy2 = matrix_mult_mm(kTe1, qx2qy2);
-  auto rhoCi = 1 / (rho * C);
   auto hvfft_kTe_qx2qy2 = matrix_sub_mm(hvfft_matrix, kTe1qx2qy2);
+  auto rhoCi = 1 / (rho * C);
   auto dTe_dt_fft = matrix_mult_sm(rhoCi, hvfft_kTe_qx2qy2);
 
   // step 3: backward fft
   auto dTe_dt = FFT::itransform(dTe_dt_fft);
 
   // step 4: temperatures_matrix+1
   auto Dt_dTedt = matrix_mult_sm(dt, dTe_dt);
   auto temperatures_matrix_1 = matrix_add_mm(temperatures_matrix, Dt_dTedt);
 
   updateTemperatures(system, temperatures_matrix_1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void ComputeTemperature::updateTemperatures(
     System &system, Matrix<complex> temperatures_matrix_1) {
   UInt size = std::sqrt(system.getNbParticles());
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < size; j++) {
       auto &part = system.getParticle(i * size + j);
       auto &matpoint = static_cast<MaterialPoint &>(part);
       auto &temperature = matpoint.getTemperature();
       temperature = temperatures_matrix_1(i, j).real();
     }
   }
 }
diff --git a/exercice_3/src/fft.hh b/exercice_3/src/fft.hh
index 0601d62..c757a8a 100644
--- a/exercice_3/src/fft.hh
+++ b/exercice_3/src/fft.hh
@@ -1,118 +1,140 @@
 #ifndef FFT_HH
 #define FFT_HH
 /* ------------------------------------------------------ */
 #include "matrix.hh"
 #include "my_types.hh"
 #include <fftw3.h>
 /* ------------------------------------------------------ */
 
 struct FFT {
 
   static Matrix<complex> transform(Matrix<complex> &m);
   static Matrix<complex> itransform(Matrix<complex> &m);
 
   static Matrix<std::complex<int>> computeFrequencies(int size);
   static Matrix<std::complex<Real>> computeFrequenciesqx2qy2(int size);
 };
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::transform(Matrix<complex> &m_in) {
   fftw_plan p;
   int N = m_in.rows(); // same as int n1 = m_in.cols();
   Matrix<complex> m_out(N);
 
   // reinterpret cast
   complex *c = m_in.data();
   auto *c_in = reinterpret_cast<fftw_complex *>(c);
 
   c = m_out.data();
   auto *c_out = reinterpret_cast<fftw_complex *>(c);
 
   p = fftw_plan_dft_2d(N, N, c_in, c_out, FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::itransform(Matrix<complex> &m_in) {
   fftw_plan p;
   int N = m_in.rows(); // same as int n1 = m_in.cols();
   Matrix<complex> m_out(N);
 
   // reinterpret cast
   complex *c = m_in.data();
   auto *c_in = reinterpret_cast<fftw_complex *>(c);
 
   c = m_out.data();
   auto *c_out = reinterpret_cast<fftw_complex *>(c);
 
   p = fftw_plan_dft_2d(N, N, c_in, c_out, FFTW_BACKWARD, FFTW_ESTIMATE);
 
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   double rescale = 1. / (N * N);
   for (auto &&entry : index(m_out)) {
     auto &val = std::get<2>(entry);
     val *= rescale;
   }
 
   return m_out;
 }
 
 /* ------------------------------------------------------ */
 
 /* ------------------------------------------------------ */
 /* implementation of computeFrequencies (not sure about the related test) */
 inline Matrix<std::complex<int>> FFT::computeFrequencies(int size) {
+  /**
+    val = 1.0 / (n * d)
+    results = empty(n, int)
+    N = (n-1)//2 + 1
+    p1 = arange(0, N, dtype=int)
+    results[:N] = p1
+    p2 = arange(-(n//2), 0, dtype=int)
+    results[N:] = p2
+    return results * val
+  */
   Matrix<std::complex<int>> frequency_matrix(size);
   double size_sample_frequencies = size / 2;
   for (auto &&entry : index(frequency_matrix)) {
-    unsigned int i = std::get<0>(entry);
-    unsigned int j = std::get<1>(entry);
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
 
-    int i_val = i;
-    if (i >= size_sample_frequencies) {
-      i_val = i_val - size;
-    }
-    int j_val = j;
-    if (j >= size_sample_frequencies) {
-      j_val = j_val - size;
-    }
-    val.real(i_val);
-    val.imag(j_val);
+    if (j < size_sample_frequencies)
+      val.real(j);
+    else
+      val.real(j - size);
+    val.imag(0);
   }
-
   return frequency_matrix;
 }
+
 /* implementation of (qx2 + qy2) (not sure about the related test) */
 
 inline Matrix<std::complex<Real>> FFT::computeFrequenciesqx2qy2(int size) {
-  Matrix<std::complex<Real>> frequency_matrix(size);
-  double size_sample_frequencies = (size) / 2;
-  for (auto &&entry : index(frequency_matrix)) {
-    unsigned int i = std::get<0>(entry);
-    unsigned int j = std::get<1>(entry);
+  /**
+     Matrix<std::complex<Real>> frequency_matrix(size);
+     double size_sample_frequencies = (size) / 2;
+     for (auto &&entry : index(frequency_matrix)) {
+     unsigned int i = std::get<0>(entry);
+     unsigned int j = std::get<1>(entry);
+     auto &val = std::get<2>(entry);
+
+     int i_val = i;
+     if (i >= size_sample_frequencies) {
+     i_val = i_val - size;
+     }
+     int j_val = j;
+     if (j >= size_sample_frequencies) {
+     j_val = j_val - size;
+     }
+     val.real(i_val * i_val + j_val * j_val);
+     val.imag(0.0);
+     }
+
+     return frequency_matrix;
+  */
+  Matrix<std::complex<int>> frequency_matrix = computeFrequencies(size);
+
+  Matrix<std::complex<Real>> qx2qy2(size);
+
+  double size_sample_frequencies = size / 2;
+  for (auto &&entry : index(qx2qy2)) {
+    int i = std::get<0>(entry);
+    int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
 
-    int i_val = i;
-    if (i >= size_sample_frequencies) {
-      i_val = i_val - size;
-    }
-    int j_val = j;
-    if (j >= size_sample_frequencies) {
-      j_val = j_val - size;
-    }
-    val.real(i_val * i_val + j_val * j_val);
-    val.imag(0.0);
-  }
+    std::complex<int> frequency = frequency_matrix(i, j);
 
-  return frequency_matrix;
+    val.real((i * i + j * j) / (size * size));
+  }
+  return qx2qy2;
 }
 
 #endif // FFT_HH
diff --git a/exercice_3/src/test_fft.cc b/exercice_3/src/test_fft.cc
index 6f2dcaa..a570b1e 100644
--- a/exercice_3/src/test_fft.cc
+++ b/exercice_3/src/test_fft.cc
@@ -1,86 +1,86 @@
 #include "fft.hh"
 #include "my_types.hh"
 #include <gtest/gtest.h>
 
 /*****************************************************************/
 TEST(FFT, transform) {
   UInt N = 512;
   Matrix<complex> m(N);
 
   Real k = 2 * M_PI / N;
   for (auto &&entry : index(m)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
     val = cos(k * i);
   }
 
   Matrix<complex> res = FFT::transform(m);
 
   for (auto &&entry : index(res)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
     if (std::abs(val) > 1e-10)
       std::cout << i << "," << j << " = " << val << std::endl;
 
     if (i == 1 && j == 0)
       ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10);
     else if (i == N - 1 && j == 0)
       ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10);
     else
       ASSERT_NEAR(std::abs(val), 0, 1e-10);
   }
 }
 /*****************************************************************/
 
 TEST(FFT, inverse_transform) {
   UInt N = 512;
   Matrix<complex> m(N);
 
   for (auto &&entry : index(m)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
     if (i == j)
       val = N;
     else
       val = 0;
   }
 
   Matrix<complex> res = FFT::itransform(m);
 
   for (auto &&entry : index(res)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
 
     if (i == 0 && j == 0)
       ASSERT_NEAR(std::abs(val), 1, 1e-10);
     else if (i > 0 && i == N - j)
       ASSERT_NEAR(std::abs(val), 1, 1e-10);
     else
       ASSERT_NEAR(std::abs(val), 0, 1e-10);
   }
 }
 /*****************************************************************/
 
 TEST(FFT, compute_frequency) {
   UInt N = 512;
   Matrix<std::complex<int>> frequency_matrix = FFT::computeFrequencies(N);
 
   for (auto &&entry : index(frequency_matrix)) {
     int i = std::get<0>(entry);
     int j = std::get<1>(entry);
     auto &val = std::get<2>(entry);
 
-    if (i == 0)
+    if (j == 0)
       ASSERT_NEAR(val.real(), 0, 1e-10);
     else if (j == 0)
       ASSERT_NEAR(val.imag(), 0, 1e-10);
-    else if (i == N - 1)
+    else if (j == N - 1)
       ASSERT_NEAR(val.real(), -1, 1e-10);
     else if (j == N - 1)
       ASSERT_NEAR(val.imag(), -1, 1e-10);
   }
 }
diff --git a/exercice_3/src/test_temperature.cc b/exercice_3/src/test_temperature.cc
index 7e8bb22..7e085cb 100644
--- a/exercice_3/src/test_temperature.cc
+++ b/exercice_3/src/test_temperature.cc
@@ -1,41 +1,119 @@
 #include "compute_temperature.hh"
 #include "material_points_factory.hh"
 #include "my_types.hh"
 #include "system_evolution.hh"
 #include <gtest/gtest.h>
 
 /*****************************************************************/
 TEST(COMPUTE_TEMPERATURE, homogeneous) {
   UInt N = 512;
 
   System system = System();
 
-  for (int i = 0; i < N; i++) {
-    for (int j = 0; j < N; j++) {
+  for (UInt i = 0; i < N; i++) {
+    for (UInt j = 0; j < N; j++) {
       // auto &val = std::get<2>(entry);
       std::stringstream sstr;
       // position (i, j, 0);  temperature 10; 0 for the rest
       std::string line =
-          std::to_string(i) + " " + std::to_string(i) + " 0 0 0 0 0 0 0 0 10 0";
+          std::to_string(i) + " " + std::to_string(j) + " 0 0 0 0 0 0 0 0 10 0";
       sstr << line;
       auto p = MaterialPointsFactory::getInstance().createParticle();
       sstr >> *p;
       system.addParticle(std::move(p));
     }
   }
 
   ComputeTemperature compute_temperature = ComputeTemperature();
   // test 5 iterations - nothing should change regardless of the number of
   // iterations
   int nsteps = 5;
   for (int step = 0; step < nsteps; ++step)
     compute_temperature.compute(system);
 
   UInt nb_particles = system.getNbParticles();
 
   for (UInt p = 0; p < nb_particles; ++p) {
     MaterialPoint *mpoint =
         static_cast<MaterialPoint *>(&system.getParticle(p));
     ASSERT_NEAR(mpoint->getTemperature(), 10, 1e-10);
   }
 }
+
+/*****************************************************************/
+TEST(COMPUTE_TEMPERATURE, volumetric_heat) {
+  // In this test, `N` has to be an int (not UInt) so that products with a
+  // negative int can yield the appropriate result
+  int N = 128;
+
+  System system = System();
+
+  int hv;
+  Real temp;
+  // Same as for `N`: `i` has to be an int (not UInt) so that products with a
+  // negative int can yield the appropriate result
+  for (int i = 0; i < N; i++) {
+    if (i * 4 < N) {
+      hv = 0;
+      temp = -2 * i;
+    } else if (i * 4 == N) {
+      hv = -N;
+      temp = -N / 2; // -2 * i;
+    } else if (i * 4 <= 3 * N) {
+      hv = 0;
+      temp = 2 * i - N;
+    } else if (i * 4 == 3 * N) {
+      hv = N;
+      temp = N / 2; // 2 * i - N;
+    } else {
+      hv = 0;
+      temp = 2 * (-i + N);
+    }
+    for (UInt j = 0; j < N; j++) {
+      // auto &val = std::get<2>(entry);
+      std::stringstream sstr;
+      // position (i, j, 0);  temperature 10; 0 for the rest
+      std::string line = std::to_string(i) + " " + std::to_string(j) + " 0" +
+                         " 0 0 0 0 0 0 0 " + std::to_string(temp) + " " +
+                         std::to_string(hv);
+      sstr << line;
+      auto p = MaterialPointsFactory::getInstance().createParticle();
+      sstr >> *p;
+      system.addParticle(std::move(p));
+    }
+  }
+
+  ComputeTemperature compute_temperature = ComputeTemperature();
+  // test 5 iterations - nothing should change regardless of the number of
+  // iterations
+  int nsteps = 100;
+  // CsvWriter csv_writer = CsvWriter("step-" + std::to_string(0) + ".csv");
+  // csv_writer.write(system);
+  for (int step = 1; step < nsteps; ++step) {
+    // csv_writer = CsvWriter("step-" + std::to_string(step) + ".csv");
+    compute_temperature.compute(system);
+    // csv_writer.write(system);
+  }
+
+  // Iterating over the systems' particles in the context of this test would be
+  // extremely inefficient since i will be the same throughout `N` consecutive
+  // iterations. We will iterate over rows and cols instead.
+
+  // Same as for `N`: `i` has to be an int (not UInt) so that
+  // products with a negative int can yield the appropriate result
+  for (int i = 0; i < N; i++) {
+    if (i * 4 <= N) {
+      temp = -2 * i;
+    } else if (i * 4 <= 3 * N) {
+      temp = 2 * i - N;
+    } else {
+      temp = 2 * (-i + N);
+    }
+    for (UInt j = 0; j < N; j++) {
+      UInt p = N * i + j;
+      MaterialPoint *mpoint =
+          static_cast<MaterialPoint *>(&system.getParticle(p));
+      ASSERT_NEAR(mpoint->getTemperature(), temp, 1);
+    }
+  }
+}