diff --git a/Correction/Makefile b/Correction/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..c6e329fd2d0139683c0bc54eac1dc87ae16ad17c
--- /dev/null
+++ b/Correction/Makefile
@@ -0,0 +1,15 @@
+.PHONY:clean all
+
+CC=nvc++
+FLAGS=-std=c++20 -stdpar=gpu -O3 -g
+
+all: artefacts/heat_gpu
+
+artefacts/%: artefacts/%.o
+	$(CC) $(FLAGS) $^ -o $@
+
+artefacts/%.o: %.cpp
+	$(CC) $(FLAGS) -c $^ -o $@
+
+clean:
+	rm -vf artefacts/**
\ No newline at end of file
diff --git a/Correction/artefacts/.gitkeep b/Correction/artefacts/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Correction/figure.py b/Correction/figure.py
new file mode 100644
index 0000000000000000000000000000000000000000..06cd53409e516f7f89c771e05afa36511e01128c
--- /dev/null
+++ b/Correction/figure.py
@@ -0,0 +1,20 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+import sys
+
+if len(sys.argv) < 3:
+    print("Usage: python script.py input_csv_file output_pdf_file")
+    sys.exit(1)
+
+file_path = sys.argv[1]
+output_file_path = sys.argv[2]
+
+data = pd.read_csv(file_path, header=None)  
+
+plt.figure(figsize=(10, 8))  
+plt.imshow(data, cmap='hot', interpolation='nearest') 
+plt.colorbar() 
+
+plt.savefig(output_file_path, bbox_inches='tight')
+
+plt.close()
\ No newline at end of file
diff --git a/Correction/heat_gpu.cpp b/Correction/heat_gpu.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea2d248af1194da7d22639149c359e1e30619e14
--- /dev/null
+++ b/Correction/heat_gpu.cpp
@@ -0,0 +1,104 @@
+#include <algorithm>
+#include <chrono>
+#include <execution>
+#include <iostream>
+#include <ranges>
+#include <string>
+
+#include "matrix.h"
+
+constexpr auto kPolicy = std::execution::par_unseq;
+
+template <typename T>
+void stepSlow(Matrix<T> &current, Matrix<T> &next) {
+    auto cols = current.cols();
+
+    auto idxs_rows = std::views::iota(1, rows - 1);
+    std::for_each(kPolicy,
+                  idxs_rows.begin(), idxs_rows.end(),
+                  [cur = current.dataSpan(),
+                   nxt = next.dataSpan(),
+                   cols](int i) {
+                      auto idxs_cols = std::views::iota(1, cols - 1);
+                      std::for_each(
+                          idxs_cols.begin(), idxs_cols.end(),
+                          [cur, nxt, cols, row = i * cols](int j) {
+                              auto idx = row + j;  // we need a 1d index for the spans
+                              nxt[idx] = 0.25 * (cur[idx - 1] +
+                                                 cur[idx + 1] +
+                                                 cur[idx - cols] +
+                                                 cur[idx + cols]);
+                          });
+                  });
+}
+
+template <typename T>
+void step(Matrix<T> &current, Matrix<T> &next) {
+    auto cols = current.cols();
+    auto rows = current.rows();
+
+    auto up = -1 * cols;
+    auto down = cols;
+    auto left = -1;
+    auto right = 1;
+
+    auto idxs = std::views::iota(0, (cols - 2) * (rows - 2));
+
+    std::for_each(kPolicy,
+                  idxs.begin(), idxs.end(),
+                  [cur = current.dataSpan(),
+                   nxt = next.dataSpan(),
+                   cols,
+                   up, down, left, right](int idx) {
+                      // We compute the 1d index to access the spans
+                      auto r = idx / (cols - 2);
+                      auto c = idx % (cols - 2);
+                      auto i = (r + 1) * cols + c + 1;
+                      nxt[i] = 0.25 * (cur[i + left] +
+                                       cur[i + right] +
+                                       cur[i + up] +
+                                       cur[i + down]);
+                  });
+}
+
+int main(int argc, char **argv) {
+    int m = std::stoi(argv[1]);
+    int n = std::stoi(argv[2]);
+    int tmax = std::stoi(argv[3]);
+    std::string filename = argv[4];
+
+    auto start_global = std::chrono::steady_clock::now();
+
+    std::cout << "processors : seq" << std::endl;
+    std::cout << "m : " << m << std::endl;
+    std::cout << "n : " << n << std::endl;
+    std::cout << "tmax : " << tmax << std::endl;
+
+    Matrix<double> *U = new Matrix<double>(m, n);
+    Matrix<double> *Unext = new Matrix<double>(m, n);
+    init_matrix(*U);
+    init_matrix(*Unext);
+
+    auto start_iter = std::chrono::steady_clock::now();
+    for (int t = 0; t < tmax; t++) {
+        step(*U, *Unext);
+        std::swap(U, Unext);
+        if (t % 100 == 0) {
+            std::cout << "t:" << t << std::endl;
+        }
+    }
+    auto end_iter = std::chrono::steady_clock::now();
+    auto duration_iter = std::chrono::duration<double>(end_iter - start_iter);
+
+    std::cout << "iterations time : " << duration_iter.count() << "[s]" << std::endl;
+
+    write_mat(*U, filename);
+    std::cout << "sum of elements of matrix : " << sum_mat(*U) << std::endl;
+
+    delete U;
+    delete Unext;
+
+    auto end_global = std::chrono::steady_clock::now();
+    auto duration_global = std::chrono::duration<double>(end_global - start_global);
+    std::cout << "total time : " << duration_global.count() << "[s]" << std::endl;
+}
\ No newline at end of file
diff --git a/Correction/matrix.h b/Correction/matrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..46fbc1c1b6d1a067970b9b898758423fcb6a8a9c
--- /dev/null
+++ b/Correction/matrix.h
@@ -0,0 +1,112 @@
+#include <fstream>
+#include <span>
+#include <stdexcept>
+
+// a simple generic Matrix class
+// templates in c++ allows genericity
+// https://en.cppreference.com/w/cpp/language/templates
+template <typename T>
+class Matrix {
+   public:
+    // constructor, allocate linear array and
+    // arrange pointer for 2D array
+    Matrix<T>(int m, int n) {
+        this->m = m;
+        this->n = n;
+        mat = new T *[m];
+        mat[0] = new T[m * n];
+        for (int i = 0; i < m; i++) {
+            mat[i] = &(mat[0][i * n]);
+        }
+    }
+    // accessor [] of class Matrix returns a pointer
+    // to the corresponding line
+    T *operator[](int index) {
+        if (index < 0 || index >= m) {
+            throw std::out_of_range("Index out of range");
+        }
+        return mat[index];
+    }
+    // destructor of class Matrix
+    // automatically called when object goes out of scope
+    ~Matrix() {
+        delete[] mat[0];
+        delete[] mat;
+    }
+    int cols() {
+        return n;
+    }
+    int rows() {
+        return m;
+    }
+
+    std::span<T> dataSpan() {
+        return std::span<T>(this->mat[0], n * m);
+    }
+
+   private:
+    int m, n;
+    T **mat;
+};
+
+// initialise a Matrix
+// here, & is the symbol for reference
+// which allows passing an object by reference
+// without using pointer syntax
+// https://en.cppreference.com/w/cpp/language/reference
+template <typename T>
+void init_matrix(Matrix<T> &mat) {
+    for (int i = 0; i < mat.rows(); i++) {
+        for (int j = 0; j < mat.cols(); j++) {
+            mat[i][j] = 0;
+        }
+    }
+    for (int i = 0; i < mat.rows(); i++) {
+        mat[i][0] = 0;
+        mat[i][mat.cols() - 1] = 1;
+    }
+    for (int j = 0; j < mat.cols(); j++) {
+        mat[0][j] = 1;
+        mat[mat.rows() - 1][j] = 0;
+    }
+}
+
+template <typename T>
+void print_mat(Matrix<T> &m) {
+    for (int i = 0; i < m.rows(); i++) {
+        for (int j = 0; j < m.cols(); j++) {
+            std::cout << m[i][j] << ", ";
+        }
+        std::cout << std::endl;
+    }
+}
+
+template <typename T>
+void write_mat(Matrix<T> &m, std::string filename) {
+    std::ofstream file(filename);
+    if (!file.is_open()) {
+        std::cerr << "Failed to open file for writing: " << filename << std::endl;
+        return;
+    }
+    for (int i = 0; i < m.rows(); i++) {
+        for (int j = 0; j < m.cols(); j++) {
+            file << m[i][j];
+            if (j != m.cols() - 1) {
+                file << ",";
+            }
+        }
+        file << std::endl;
+    }
+    file.close();
+}
+
+template <typename T>
+T sum_mat(Matrix<T> &m) {
+    T sum = 0;
+    for (int i = 0; i < m.rows(); i++) {
+        for (int j = 0; j < m.cols(); j++) {
+            sum += m[i][j];
+        }
+    }
+    return sum;
+}
\ No newline at end of file
diff --git a/Correction/run_baobab.sh b/Correction/run_baobab.sh
new file mode 100644
index 0000000000000000000000000000000000000000..3dd28fba7bc1c8677b1f877e342d1589f5ad8101
--- /dev/null
+++ b/Correction/run_baobab.sh
@@ -0,0 +1,11 @@
+#!/bin/sh
+#SBATCH --job-name=heat_gpu
+#SBATCH --output=heat_gpu.o%j
+#SBATCH --time=00:30:00
+#SBATCH --partition=shared-gpu
+#SBATCH --gpus=1
+#SBATCH --constraint="COMPUTE_TYPE_AMPERE|COMPUTE_TYPE_TURING|"
+
+module load NVHPC/24.1-CUDA-12.3.0
+
+srun ./artefacts/heat_gpu 10000 10000 10000 output_$SLURM_JOBID.csv
\ No newline at end of file
diff --git a/Correction/run_one_job.sh b/Correction/run_one_job.sh
new file mode 100644
index 0000000000000000000000000000000000000000..c26c68898eb76ccf6aa10240cd42ff0685faa275
--- /dev/null
+++ b/Correction/run_one_job.sh
@@ -0,0 +1,16 @@
+#!/bin/sh
+#SBATCH --job-name=heat
+#SBATCH --output=heat.o%j
+#SBATCH --cpus-per-task=1
+#SBATCH --mem-per-cpu=2000
+#SBATCH --partition=shared-cpu
+#SBATCH --time=01:00:00
+#SBATCH --constraint=EPYC-7742
+
+echo $SLURM_NODELIST
+
+module purge
+module load foss/2023a
+
+srun heat_par 10000 10000 10000 output_$SLURM_JOBID.csv
+rm output_$SLURM_JOBID.csv
diff --git a/Correction/run_yggdrasil.sh b/Correction/run_yggdrasil.sh
new file mode 100644
index 0000000000000000000000000000000000000000..5fef0d22cf1d18183bac5e51232aa7bc69ae152a
--- /dev/null
+++ b/Correction/run_yggdrasil.sh
@@ -0,0 +1,11 @@
+#!/bin/sh
+#SBATCH --job-name=heat_gpu
+#SBATCH --output=heat_gpu.o%j
+#SBATCH --time=00:30:00
+#SBATCH --partition=shared-gpu
+#SBATCH --gpus=1
+#SBATCH --constraint="COMPUTE_TYPE_TURING"
+
+module load NVHPC/24.1-CUDA-12.3.0
+
+srun ./artefacts/heat_gpu 10000 10000 10000 output_$SLURM_JOBID.csv
\ No newline at end of file