diff --git a/src/numpy_pandas/numerical_methods.py b/src/numpy_pandas/numerical_methods.py
index ff90bb9..ed72637 100644
--- a/src/numpy_pandas/numerical_methods.py
+++ b/src/numpy_pandas/numerical_methods.py
@@ -93,23 +93,33 @@ def sieve_of_eratosthenes(n: int) -> List[int]:
 
 def linear_equation_solver(A: List[List[float]], b: List[float]) -> List[float]:
     n = len(A)
-    augmented = [row[:] + [b[i]] for i, row in enumerate(A)]
+    augmented = [A[i][:] + [b[i]] for i in range(n)]
+    # Gaussian Elimination with Partial Pivoting
     for i in range(n):
-        max_idx = i
+        # Find pivot
+        max_idx, max_value = i, abs(augmented[i][i])
         for j in range(i + 1, n):
-            if abs(augmented[j][i]) > abs(augmented[max_idx][i]):
-                max_idx = j
-        augmented[i], augmented[max_idx] = augmented[max_idx], augmented[i]
+            val = abs(augmented[j][i])
+            if val > max_value:
+                max_idx, max_value = j, val
+        if max_idx != i:
+            augmented[i], augmented[max_idx] = augmented[max_idx], augmented[i]
+        ai = augmented[i]
+        inv_aii = 1.0 / ai[i]
         for j in range(i + 1, n):
-            factor = augmented[j][i] / augmented[i][i]
+            rowj = augmented[j]
+            factor = rowj[i] * inv_aii
+            # In-place update
             for k in range(i, n + 1):
-                augmented[j][k] -= factor * augmented[i][k]
-    x = [0] * n
+                rowj[k] -= factor * ai[k]
+    # Back substitution
+    x = [0.0] * n
     for i in range(n - 1, -1, -1):
-        x[i] = augmented[i][n]
+        ai = augmented[i]
+        sum_ax = 0.0
         for j in range(i + 1, n):
-            x[i] -= augmented[i][j] * x[j]
-        x[i] /= augmented[i][i]
+            sum_ax += ai[j] * x[j]
+        x[i] = (ai[n] - sum_ax) / ai[i]
     return x