diff --git a/src/language/algorithms/ParabolicHeat.cpp b/src/language/algorithms/ParabolicHeat.cpp
index 5dcfb1bfeb186fab3693ffb32697f755c067aa4b..4a0857bfa858e002e1131e381cb725dcfe973174 100644
--- a/src/language/algorithms/ParabolicHeat.cpp
+++ b/src/language/algorithms/ParabolicHeat.cpp
@@ -235,8 +235,9 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
 
     CellValue<double> Tj =
       InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-    CellValue<double> Temperature{mesh->connectivity()};
-    Temperature.fill(1);
+    CellValue<double> Temperature =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
+    //{mesh->connectivity()};
     FaceValue<double> Tl =
       InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
     FaceValue<double> Temperature_face =
@@ -254,7 +255,14 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
     }
     double time = 0;
     Assert(dt > 0, "The time-step have to be positive!");
+    const CellValue<const double> primal_Vj = mesh_data.Vj();
     do {
+      {
+        VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension) + std::to_string(time), 0.01);
+
+        vtk_writer.write(mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj}}, 0,
+                         true);   // forces last output
+      }
       double deltat = std::min(dt, Tf - time);
       std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
       ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face, Tf,
@@ -265,10 +273,9 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
       Vector<double> error{mesh->numberOfCells()};
       CellValue<double> cell_error{mesh->connectivity()};
       Vector<double> face_error{mesh->numberOfFaces()};
-      double error_max                        = 0.;
-      size_t cell_max                         = 0;
-      const CellValue<const double> primal_Vj = mesh_data.Vj();
-      const FaceValue<const double> mes_l     = [&] {
+      double error_max                    = 0.;
+      size_t cell_max                     = 0;
+      const FaceValue<const double> mes_l = [&] {
         if constexpr (Dimension == 1) {
           FaceValue<double> compute_mes_l{mesh->connectivity()};
           compute_mes_l.fill(1);