diff --git a/src/nabladelta.c b/src/nabladelta.c
new file mode 100644
index 0000000..0c69e65
--- /dev/null
+++ b/src/nabladelta.c
@@ -0,0 +1,100 @@
+#include "nabladelta.h"
+#include "matrix_operation.h"
+#include <stdlib.h>
+#include <math.h>
+#include <float.h>
+#include <stdio.h>
+
+double* nabladelta(const SparseMatrix *matrix, double epsilon) {
+    int N = matrix->num_nodes;
+    int num_arcs = matrix->num_arcs;
+
+    double *nabla = malloc(N * sizeof(double));
+    double *delta = malloc(N * sizeof(double));
+    double *X = malloc(N * sizeof(double));
+    double *Y = malloc(N * sizeof(double));
+    double *temp = malloc(N * sizeof(double));
+
+    if (!nabla || !delta || !X || !Y || !temp) {
+        fprintf(stderr, "Memory allocation failed\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // Initialize nabla, delta
+    // ∇[j] = min(i, G[i, j]),  nabla(j) = minimum of column j
+    // ∆[j] = max(i, G[i, j]),  delta(j) = maximum of column j
+
+    for (int j = 0; j < N; j++) {
+        nabla[j] = 0.0;
+        delta[j] = 0.0;
+    }
+
+    for (int i = 0; i < num_arcs; i++) {
+        int j = matrix->arcs[i].dest;
+        double value = matrix->arcs[i].value;
+
+        if (value < nabla[j]) {
+            nabla[j] = value;
+        }
+        if (value > delta[j]) {
+            delta[j] = value;
+        }
+    }
+
+    // 1. X(0) ← ∇
+    // 2. Y (0) ← ∆
+    for (int i = 0; i < N; i++) {
+        X[i] = nabla[i];
+        Y[i] = delta[i];
+    }
+
+    // 3.
+    double diff = 1.0;
+    int iter = 0;
+    while (diff > epsilon) {
+        // ||X|| is norm 1
+        double X_norm1 = 0.0;
+        double Y_norm1 = 0.0;
+        for (int i = 0; i < N; i++) {
+            X_norm1 += X[i];
+            Y_norm1 += Y[i];
+        }
+        printf("Iteration %d: ||X|| = %f\n", iter, X_norm1);
+        printf("Iteration %d: ||Y|| = %f\n", iter, Y_norm1);
+
+
+        // (a): X(k+1) = max(X(k), X(k)G + nabla(1 - ||X(k)||))
+        multiply_vector_matrix_parallel(X, matrix, temp);  // X(k)*G
+
+        for (int i = 0; i < N; i++) {
+            double update = temp[i] + nabla[i] * (1.0 - X_norm1);
+            X[i] = fmax(X[i], update);
+        }
+
+        // (b): Y(k+1) = min(Y(k), Y(k)G + Delta(1 - ||Y(k)||))
+        multiply_vector_matrix_parallel(Y, matrix, temp);  // Y(k)*G
+
+        for (int i = 0; i < N; i++) {
+            double update = temp[i] + delta[i] * (1.0 - Y_norm1);
+            Y[i] = fmin(Y[i], update);
+        }
+
+        // (c): ||X(k) - Y(k)|| < ε
+        diff = 0.0;
+        for (int i = 0; i < N; i++) {
+            diff += fabs(X[i] - Y[i]);
+        }
+
+        if ((++iter)%5 == 0) {
+            printf("Iteration %d: diff = %.16f\n", iter, diff);
+        }
+    }
+
+    free(nabla);
+    free(delta);
+    free(Y);
+    free(temp);
+
+    // X~pi~Y (convergence)
+    return X;
+}
diff --git a/src/nabladelta.h b/src/nabladelta.h
new file mode 100644
index 0000000..6163306
--- /dev/null
+++ b/src/nabladelta.h
@@ -0,0 +1,8 @@
+#ifndef NABLADELTA_H
+#define NABLADELTA_H
+
+#include "sparse_matrix.h"
+
+double* nabladelta(const SparseMatrix *matrix, double epsilon);
+
+#endif