diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp
index a0cf9d6..2c73633 100644
--- a/mgrid_2/getmatrix.cpp
+++ b/mgrid_2/getmatrix.cpp
@@ -394,6 +394,9 @@ void FEM_Matrix::CalculateLaplace_mult(vector<double> &f)
     for (int k = 0; k < _nrows; ++k) {
         _sk[k] = 0.0;
     }
+    for (int k = 0; k < _nrows; ++k) {
+        f[k] = 0.0;
+    }
 
     double ske[3][3], fe[3];
     //  Loop over all elements
@@ -457,17 +460,17 @@ double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
     switch (subdomain)
     {
         // ceramic mug
-        case 1:
+        case 0:
             c = 2.0 * 1e6;
             break;
 
         // water
-        case 2:
+        case 1:
             c = 4.184 * 1e6;       
             break;
 
         // air
-        case 3:
+        case 2:
             c = 1.2 * 1e3;
             break;
 
@@ -479,7 +482,7 @@ double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
     return c;
 }
 
-void FEM_Matrix::AddMass_mult(vector<double> &f)
+void FEM_Matrix::AddMass_mult(vector<double> &f, const double scale_factor)
 {
     cout << "\n############   FEM_Matrix::AddMass_mult ";
     double tstart = omp_get_wtime();                  // OpenMP
@@ -503,7 +506,7 @@ void FEM_Matrix::AddMass_mult(vector<double> &f)
 
         //cout << subdomain << endl;
 
-        CalcElem_MasseSpecific(ia.data() + 3 * i, xc.data(), c, ske);
+        CalcElem_MasseSpecific(ia.data() + 3 * i, xc.data(), c * scale_factor, ske);
         //AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
         AddElem_3(ia.data() + 3 * i, ske, fe, f);
     }
diff --git a/mgrid_2/getmatrix.h b/mgrid_2/getmatrix.h
index ac7dcc9..1f410be 100644
--- a/mgrid_2/getmatrix.h
+++ b/mgrid_2/getmatrix.h
@@ -365,7 +365,7 @@ class FEM_Matrix: public CRS_Matrix
         *
         * @param[in,out] f (preallocated) rhs/load vector
         */
-      void AddMass_mult(std::vector<double> &f);
+      void AddMass_mult(std::vector<double> &f, const double scale_factor);
 
         /**
         * Calculates the entries of f.e. stiffness matrix for the Laplace operator
diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp
index 7dddf2a..91eb1b3 100644
--- a/mgrid_2/main.cpp
+++ b/mgrid_2/main.cpp
@@ -26,6 +26,18 @@ int main(int argc, char **argv )
     //mesh.Debug();
     //mesh_c.DebugEdgeBased();
 
+    // ##########################################
+    // Parameteres
+    // ##########################################
+
+    double dt = 1.0;            // time step
+    int steps = 1;            // number of time iterations
+
+    double u0_mug = 18.0;
+    double u0_fluid = 80.0;
+    double u0_air = 18.0;
+    double u_out = 18.0;
+
     // ##########################################
     // Assembling
     // ##########################################
@@ -33,40 +45,53 @@ int main(int argc, char **argv )
     // Initializing FEM matrix !pattern! (only zero entries now)
     FEM_Matrix SK(mesh_c);                                    // CRS matrix
     //SK.writeBinary("sparseMatrix.bin");
-    //SK.Debug();
+    // SK.Debug();
 
     vector<double> fv(SK.Nrows(), 0.0);                    
-    SK.CalculateRHS(fv, [](double x, double y) {return 0;});  // r.h.s.
+    SK.CalculateRHS(fv, [](double x, double y) {return 0;});  // rhs (f=0)
 
-    SK.CalculateLaplace_mult(fv);                             // stiffness matrix
-    SK.AddMass_mult(fv);                                      // mass matrix
-    SK.ApplyRobinBC_mult(fv, 18.0);                           // apply Robin bnd
+    SK.CalculateLaplace_mult(fv);                             // stiffness matrix (+K)
+    SK.AddMass_mult(fv, 1.0/dt);                              // add mass matrix (+M/dt)
+    SK.ApplyRobinBC_mult(fv, u_out);                           // apply Robin-bnd (+C = +F)
+                                                              // SK = M/dt + K + C = F
+    
+    // SK.Debug();
+    // SK.CheckRowSum();
+    // SK.CheckMatrix();
+    
+    FEM_Matrix Mdt(mesh_c);
+    Mdt.AddMass_mult(fv, 1.0/dt);                              // mass matrix (Mdt = M/dt)
 
-    SK.CheckRowSum();
-    SK.CheckMatrix();
+    // Mdt.Debug();
+    // Mdt.CheckRowSum();
+    // Mdt.CheckMatrix();
 
     // ##########################################
-    // Timestepping
+    // Timestepping        (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
     // ##########################################
 
-    double dt = 1.0;            // time step
-    int steps = 100;            // number of time iterations
-
-    // Initialize temperature
-    vector<double> uv(SK.Nrows(), 0.0);    // temperature
-    mesh_c.Init_Solution_mult(uv, 0, [](double x, double y) -> double { return 18; });    // mug
-    mesh_c.Init_Solution_mult(uv, 1, [](double x, double y) -> double { return 80; });    // fluid
-    mesh_c.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; });    // air
-
-    for (int i = 0; i < steps; ++i)
-    {
-        // TODO        
-    }
+    // Initialize temperature u_0
+    vector<double> uv(SK.Nrows(), 0.0);                                                               // temperature
+    mesh_c.Init_Solution_mult(uv, 0, [u0_mug](double x, double y) -> double { return u0_mug; });      // mug
+    mesh_c.Init_Solution_mult(uv, 1, [u0_fluid](double x, double y) -> double { return u0_fluid; });  // fluid
+    mesh_c.Init_Solution_mult(uv, 2, [u0_air](double x, double y) -> double { return u0_air; });      // air
 
+    // TODO DINO
     auto t3 = system_clock::now(); // start timer
-    //JacobiSolve(SK, fv, uv );          // solve the system of equations
+    for (int step = 0; step < steps; ++step)
+    {
+        vector<double> G(Mdt.Nrows(), 0.0);        
+        Mdt.Mult(G, uv);                                           // G = M/dt * u_{n}
+        for (size_t i = 0; i < Mdt.Nrows(); ++i)
+        {
+            fv[i] += G[i];                                         // F + G
+        }
+
+        JacobiSolve(SK, fv, uv);          // solve: (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
+    }
     auto t4 = system_clock::now();  // stop timer
 
+
     auto duration = duration_cast<microseconds>(t4 - t3);        // duration in microseconds
     double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
     cout << "JacobiSolve: timing in sec. : " << t_diff << endl;