diff --git a/Homework3/src/compute_temperature.cc b/Homework3/src/compute_temperature.cc
index b7a5fae..f5381fd 100644
--- a/Homework3/src/compute_temperature.cc
+++ b/Homework3/src/compute_temperature.cc
@@ -1,66 +1,68 @@
 #include "compute_temperature.hh"
 #include "fft.hh"
 #include "material_point.hh"
 #include <cmath>
 
 ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {}
 
 /* -------------------------------------------------------------------------- */
 
 void ComputeTemperature::setDeltaT(Real dt) {
 
     this->dt = dt;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void ComputeTemperature::compute(System& system) {
 
     UInt Nb = system.getNbParticles();
     UInt dim = UInt(sqrt(Nb));
 
     Matrix<complex> theta(dim);
     Matrix<complex> h_v(dim);
     Matrix<complex> theta_fourier(dim);
     Matrix<complex> h_v_fourier(dim);
     Matrix<std::complex<int>> wave_nb(dim);
     Matrix<complex> dthetadt(dim);
 
     // Initialize matrices
     for (auto&& par : index(theta)){
         int i = std::get<0>(par);
         int j = std::get<1>(par);
-        UInt id = i + j*dim;
-        //UInt id =  i*dim + j;
+        //UInt id = i + j*dim;
+        UInt id =  i*dim + j;
         auto& ic = dynamic_cast<MaterialPoint&>(system.getParticle(id));
-        theta(i,j) = ic.getTemperature();
-        h_v(i,j) = ic.getHeatRate();
+        theta(i,j) = complex(ic.getTemperature(),0);
+        h_v(i,j)   = complex(ic.getHeatRate()   ,0);
     }
 
     theta_fourier = FFT::transform(theta);
     h_v_fourier = FFT::transform(h_v);
     wave_nb = FFT::computeFrequencies(dim);
 
-    Real factor = 4.0*M_PI*M_PI/(L*L);
+    Real factor = pow(2.0*M_PI/L,2);
     for(int j=0; j<dim; j++){
         for(int i=0; i<dim; i++){
-            dthetadt(i, j) = h_v_fourier(i, j) - k*theta_fourier(i, j)*factor*Real(norm(wave_nb(i, j)));
-            dthetadt(i, j) = dthetadt(i, j)/(rho*C);
+            dthetadt(i, j) = -k*theta_fourier(i,j)*factor*norm(wave_nb(i,j));
+                             +h_v_fourier(i,j)
+            dthetadt(i,j) = dthetadt(i,j)/(rho*C);
         }
     }
 
     dthetadt = FFT::itransform(dthetadt);
 
     for (auto&& par : index(dthetadt)){
         int i = std::get<0>(par);
         int j = std::get<1>(par);
-        UInt id = i + j*dim;
+        //UInt id = i + j*dim;
+        UInt id =  i*dim + j;
         auto& value = std::get<2>(par);
         dynamic_cast<MaterialPoint&>(system.getParticle(id)).getTemperature() += value.real()*this->dt;
     }
 
 
 }
 
 /* -------------------------------------------------------------------------- */