Error in Back Substitution Function For Gaussian Elimination of a matrix

I'm trying to build a program to perform Gaussian Elimination on an matrix, I was able to create a function to convert the matrix into row echelon form. But my Back Substitution program has an error which I tried to solve in many ways but couldn't.

Here's my code:

def back_substitution(M):
    """
    Perform back substitution on an augmented matrix (with unique solution) in reduced row echelon form to find the solution to the linear system.

    Parameters:
    - M (numpy.array): The augmented matrix in row echelon form with unitary pivots (n x n+1).

    Returns:
    numpy.array: The solution vector of the linear system.
    """
    
    # Make a copy of the input matrix to avoid modifying the original
    M = M.copy()

    # Get the number of rows (and columns) in the matrix of coefficients
    num_rows = M.shape[0]

    if num_rows == 0:
        return np.array([])
    
    # Iterate from bottom to top
    for row in reversed(range(num_rows)): 
        
        substitution_row = M[row, :-1]  # Extract the row without the last element

        # Get the index of the first non-zero element in the substitution row
        index = np.argmax(substitution_row)

        # Iterate over the rows above the substitution_row
        for j in range(row): 

            row_to_reduce = M[j, :-1]

            # Get the value of the element at the found index in the row to reduce
            value = row_to_reduce[index]
            
            # Perform the back substitution step using the formula row_to_reduce -> row_to_reduce - value * substitution_row
            row_to_reduce -= value * substitution_row

            # Replace the updated row in the matrix
            M[j, :-1] = row_to_reduce
            
    # Extract the solution from the last column
    solution = M[:, -1]
    
    return solution

The Error:

Wrong output for test case check_matrix_1. 
    Expected:
     [-0.33333333 -1.33333333  2.33333333].
    Got:
 [4.         1.7        2.33333333].
Wrong output for test case check_matrix_2. 
    Expected:
     [0.7857143  1.64285714 0.        ].
    Got:
 [9.         1.64285714 0.        ].
Wrong output for test case check_matrix_3. 
    Expected:
     [19. -2.  1.].
    Got:
 [9. 6. 1.].
 1  Tests passed
 3  Tests failed

Thanks in advance!

I tried modifying the row_to_reduce variable in everyway I could think of but couldn't find a soltuion.

Solution

Assuming matrix M is full rank and in a row echelon form. Then you can use the following:

import numpy as np

def back_solve(M):
  n = M.shape[0]
  solution = np.zeros(n)
  for i in reversed(range(n)):
    solution[i] =  (M[i,-1] - solution @ M[i, :-1]) / M[i,i]
  return solution

EXAMPLE:

a = np.array([[1,2,3],[0,4,5],[0,0,6]])

np.linalg.solve(a, [1,2,3])
array([-0.25 , -0.125,  0.5  ])

back_solve(np.c_[a, [1,2,3]])
array([-0.25 , -0.125,  0.5  ])