diff --git a/CMakeLists.txt b/CMakeLists.txt
index 599ed481..cbbcddc1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,7 +5,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 set(DIFF_ENGINE_VERSION_MAJOR 0)
 set(DIFF_ENGINE_VERSION_MINOR 1)
-set(DIFF_ENGINE_VERSION_PATCH 4)
+set(DIFF_ENGINE_VERSION_PATCH 5)
 set(DIFF_ENGINE_VERSION "${DIFF_ENGINE_VERSION_MAJOR}.${DIFF_ENGINE_VERSION_MINOR}.${DIFF_ENGINE_VERSION_PATCH}")
 add_compile_definitions(DIFF_ENGINE_VERSION="${DIFF_ENGINE_VERSION}")
 
diff --git a/include/elementwise_univariate.h b/include/elementwise_univariate.h
index 54ade568..69cc5ab6 100644
--- a/include/elementwise_univariate.h
+++ b/include/elementwise_univariate.h
@@ -37,6 +37,7 @@ expr *new_atanh(expr *child);
 expr *new_logistic(expr *child);
 expr *new_power(expr *child, double p);
 expr *new_xexp(expr *child);
+expr *new_normal_cdf(expr *child);
 
 /* the jacobian and wsum_hess for elementwise univariate atoms are always
    initialized in the same way and implement the chain rule in the same way */
diff --git a/src/affine/broadcast.c b/src/affine/broadcast.c
index c16fe108..4d8277ea 100644
--- a/src/affine/broadcast.c
+++ b/src/affine/broadcast.c
@@ -129,9 +129,9 @@ static void jacobian_init(expr *node)
         int offset = 0;
         for (int i = 0; i < node->d2; i++)
         {
-            int nnz_in_row = Jx->p[i + 1] - Jx->p[i];
             for (int j = 0; j < node->d1; j++)
             {
+                int nnz_in_row = Jx->p[j + 1] - Jx->p[j];
                 J->p[i * node->d1 + j] = offset;
                 offset += nnz_in_row;
             }
diff --git a/src/elementwise_univariate/normal_cdf.c b/src/elementwise_univariate/normal_cdf.c
new file mode 100644
index 00000000..727a36a9
--- /dev/null
+++ b/src/elementwise_univariate/normal_cdf.c
@@ -0,0 +1,70 @@
+/*
+ * Copyright 2026 Daniel Cederberg and William Zhang
+ *
+ * This file is part of the DNLP-differentiation-engine project.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "elementwise_univariate.h"
+#include <math.h>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#ifndef M_SQRT2
+#define M_SQRT2 1.41421356237309504880
+#endif
+
+static const double INV_SQRT_2PI = 0.3989422804014326779399461;
+
+static void forward(expr *node, const double *u)
+{
+    node->left->forward(node->left, u);
+
+    double *x = node->left->value;
+    for (int i = 0; i < node->size; i++)
+    {
+        node->value[i] = 0.5 * (1.0 + erf(x[i] / M_SQRT2));
+    }
+}
+
+static void local_jacobian(expr *node, double *vals)
+{
+    double *x = node->left->value;
+    for (int j = 0; j < node->size; j++)
+    {
+        vals[j] = INV_SQRT_2PI * exp(-0.5 * x[j] * x[j]);
+    }
+}
+
+static void local_wsum_hess(expr *node, double *out, const double *w)
+{
+    double *x = node->left->value;
+    for (int j = 0; j < node->size; j++)
+    {
+        /* could avoid recomputing this (like in logistic) */
+        double phi = INV_SQRT_2PI * exp(-0.5 * x[j] * x[j]);
+        out[j] = w[j] * (-x[j] * phi);
+    }
+}
+
+expr *new_normal_cdf(expr *child)
+{
+    expr *node = new_elementwise(child);
+    node->forward = forward;
+    node->local_jacobian = local_jacobian;
+    node->local_wsum_hess = local_wsum_hess;
+
+    return node;
+}
diff --git a/tests/all_tests.c b/tests/all_tests.c
index f9fe909c..8f839f49 100644
--- a/tests/all_tests.c
+++ b/tests/all_tests.c
@@ -15,6 +15,7 @@
 #include "forward_pass/composite/test_composite.h"
 #include "forward_pass/elementwise/test_exp.h"
 #include "forward_pass/elementwise/test_log.h"
+#include "forward_pass/elementwise/test_normal_cdf.h"
 #include "forward_pass/test_left_matmul_dense.h"
 #include "forward_pass/test_matmul.h"
 #include "forward_pass/test_prod_axis_one.h"
@@ -101,6 +102,7 @@ int main(void)
     mu_run_test(test_promote_scalar_to_vector, tests_run);
     mu_run_test(test_exp, tests_run);
     mu_run_test(test_log, tests_run);
+    mu_run_test(test_normal_cdf, tests_run);
     mu_run_test(test_composite, tests_run);
     mu_run_test(test_sum_axis_neg1, tests_run);
     mu_run_test(test_sum_axis_0, tests_run);
diff --git a/tests/forward_pass/elementwise/test_normal_cdf.h b/tests/forward_pass/elementwise/test_normal_cdf.h
new file mode 100644
index 00000000..b9cc254f
--- /dev/null
+++ b/tests/forward_pass/elementwise/test_normal_cdf.h
@@ -0,0 +1,26 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "affine.h"
+#include "elementwise_univariate.h"
+#include "expr.h"
+#include "minunit.h"
+#include "test_helpers.h"
+
+#ifndef M_SQRT2
+#define M_SQRT2 1.41421356237309504880
+#endif
+
+const char *test_normal_cdf(void)
+{
+    double u[3] = {1.0, 2.0, 3.0};
+    expr *var = new_variable(3, 1, 0, 3);
+    expr *node = new_normal_cdf(var);
+    node->forward(node, u);
+    /* computed in python */
+    double correct[3] = {0.8413447460685429, 0.9772498680518208, 0.9986501019683699};
+    mu_assert("fail", cmp_double_array(node->value, correct, 3));
+    free_expr(node);
+    return 0;
+}