diff --git a/include/data_structures/Binary_search.h b/include/data_structures/Binary_search.h
new file mode 100644
index 00000000..964fd02d
--- /dev/null
+++ b/include/data_structures/Binary_search.h
@@ -0,0 +1,92 @@
+/*
+ * Copyright 2025 Daniel Cederberg
+ *
+ * This file is part of the PSLP project (LP Presolver).
+ *
+ * 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.
+ */
+
+#ifndef CORE_BINARY_SEARCH_H
+#define CORE_BINARY_SEARCH_H
+
+#define BSEARCH_LINEAR_THRESHOLD 8
+
+/* Find the index of 'target' in a sorted array 'arr' of length 'len'.
+   Returns the relative index (0-based) if found, or -1 if not found.
+   Falls back to linear scan for short arrays. */
+static inline int sorted_find(const int *arr, int len, int target)
+{
+    if (len <= BSEARCH_LINEAR_THRESHOLD)
+    {
+        for (int i = 0; i < len; ++i)
+        {
+            if (arr[i] == target) return i;
+            if (arr[i] > target) return -1;
+        }
+        return -1;
+    }
+
+    int lo = 0, hi = len - 1;
+    while (lo <= hi)
+    {
+        int mid = lo + (hi - lo) / 2;
+        if (arr[mid] == target)
+        {
+            return mid;
+        }
+        if (arr[mid] < target)
+        {
+            lo = mid + 1;
+        }
+        else
+        {
+            hi = mid - 1;
+        }
+    }
+    return -1;
+}
+
+/* Find the first index where arr[i] >= target in a sorted array of length 'len'.
+   Returns 'len' if all elements are less than target.
+   Falls back to linear scan for short arrays. */
+static inline int sorted_lower_bound(const int *arr, int len, int target)
+{
+    if (len <= BSEARCH_LINEAR_THRESHOLD)
+    {
+        for (int i = 0; i < len; ++i)
+        {
+            if (arr[i] >= target)
+            {
+                return i;
+            }
+        }
+        return len;
+    }
+
+    int lo = 0, hi = len;
+    while (lo < hi)
+    {
+        int mid = lo + (hi - lo) / 2;
+        if (arr[mid] < target)
+        {
+            lo = mid + 1;
+        }
+        else
+        {
+            hi = mid;
+        }
+    }
+    return lo;
+}
+
+#endif // CORE_BINARY_SEARCH_H
diff --git a/src/core/Matrix.c b/src/core/Matrix.c
index 81bd1568..f3a3d214 100644
--- a/src/core/Matrix.c
+++ b/src/core/Matrix.c
@@ -17,6 +17,7 @@
  */
 
 #include "Matrix.h"
+#include "Binary_search.h"
 #include "Debugger.h"
 #include "Memory_wrapper.h"
 #include "Numerics.h"
@@ -416,23 +417,26 @@ void print_row_starts(const RowRange *row_ranges, size_t len)
 
 double insert_or_update_coeff(Matrix *A, int row, int col, double val, int *row_size)
 {
-    int i, start, end, insertion;
+    // int i, start, end, insertion;
     double old_val = 0.0;
-    start = A->p[row].start;
-    end = A->p[row].end;
-    insertion = end;
+    int start = A->p[row].start;
+    int end = A->p[row].end;
+    // int insertion = end;
 
     // -----------------------------------------------------------------
     //             find where it should be inserted
     // -----------------------------------------------------------------
-    for (i = start; i < end; ++i)
-    {
-        if (A->i[i] >= col)
-        {
-            insertion = i;
-            break;
-        }
-    }
+    // for (i = start; i < end; ++i)
+    // {
+    //     if (A->i[i] >= col)
+    //     {
+    //         insertion = i;
+    //         break;
+    //     }
+    // }
+
+    int rel_ins = sorted_lower_bound(A->i + start, end - start, col);
+    int insertion = start + rel_ins;
 
     // -----------------------------------------------------------------
     // Insert the new value if it is nonzero. If it exists or should be
diff --git a/src/explorers/DtonsEq.c b/src/explorers/DtonsEq.c
index 27cf973a..adbe8292 100644
--- a/src/explorers/DtonsEq.c
+++ b/src/explorers/DtonsEq.c
@@ -18,6 +18,7 @@
 
 #include "DTonsEq.h"
 #include "Activity.h"
+#include "Binary_search.h"
 #include "Bounds.h"
 #include "Constraints.h"
 #include "CoreTransformations.h"
@@ -240,7 +241,7 @@ static inline
                                         double aij, double aik, int *row_size,
                                         PostsolveInfo *postsolve_info)
 {
-    int ii, start, end, insertion;
+    int start, end, insertion;
     double old_val = 0.0;
     int subst_idx = -1;
     int diff_row_size = 0;
@@ -252,25 +253,13 @@ static inline
     // find index of variable that is substituted and and the index of
     // insertion of the variable that stays
     // -----------------------------------------------------------------
-    for (ii = start; ii < end; ++ii)
-    {
-        if (A->i[ii] == k)
-        {
-            subst_idx = ii;
-            break;
-        }
-    }
-
-    assert(subst_idx != -1);
+    int row_len = end - start;
+    int rel = sorted_find(A->i + start, row_len, k);
+    assert(rel != -1);
+    subst_idx = start + rel;
 
-    for (ii = start; ii < end; ++ii)
-    {
-        if (A->i[ii] >= j)
-        {
-            insertion = ii;
-            break;
-        }
-    }
+    int rel_ins = sorted_lower_bound(A->i + start, row_len, j);
+    insertion = start + rel_ins;
 
     // -----------------------------------------------------------------
     // Compute new coefficient of the variable that stays.