def bwd_subs( U, y ):
    """Solve the linear system Ux = y with upper triangular matrix U by forward substitution."""
    n = U.shape[0]
    result = np.zeros(n) # initiate an array of the same length as y filled with zeros as floats
    n = n - 1 # Because Indices start at 0
    result[n] = y[n] / U[n, n] # get the last coefficient
    for j in range(n, -1, -1): # iterate backwards through the rows
        help = 0.0
        for k in range(j + 1, n + 1): # iterate through all of the known values
            help = U[j, k] * result[k] + help # sum up all of the known values
        result[j] = (1 / U[j, j]) * (y[j] - help) # formula from the script

    return result # return an array with all of the coefficients

def gauss_solve( A, b ):
    """Solve the linear system Ax=b using direct Gaussian elimination and backward substitution."""

    n = A.shape[0]
    bCopy = np.copy(b)
    n = n - 1 # Because Indices start at 0
    for i in range(0, n): # iterate through your matrix
        if A[i, i] == 0.: # if the value at this spot is 0 you cannot apply the algorithm
            for i2 in range(i + 1, n + 1): # check the following rows because you cannot pick one from before
                if A[i2, i] != 0.: # if 1 is found
                    A[[i, i2]] = A[[i2, i]] # swap the rows
                    bCopy[[i, i2]] = bCopy[[i2, i]] # swap the rows
                    break # no need to search for more
        for i3 in range(i + 1, n + 1): # perform the algorithm
            factor = A[i3, i] / A[i, i] # the times you have to subtract the first row
            A[i3] = A[i3] - (factor * A[i]) # subtract the first row factor times
            bCopy[i3] = bCopy[i3] - (factor * bCopy[i]) # same for the results

    return bwd_subs(A, bCopy) # now that they are in the right form we can solve them