Résoudre par le pivot, et détecter les trois cas

def resoudre(A, Y, tol=1e-12):
    n = len(A)
    # matrice augmentee, copiee pour ne pas modifier l'appelant
    M = [A[i][:] + [Y[i]] for i in range(n)]
    for c in range(n):
        # pivot partiel : le plus grand coefficient disponible
        p = max(range(c, n), key=lambda i: abs(M[i][c]))
        if abs(M[p][c]) < tol:
            return None                 # colonne sans pivot : pas unique
        M[c], M[p] = M[p], M[c]
        for i in range(n):
            if i != c:
                f = M[i][c] / M[c][c]
                for j in range(c, n + 1):
                    M[i][j] = M[i][j] - f * M[c][j]
    return [M[i][n] / M[i][i] for i in range(n)]

Les trois cas, sur le système $x+y+z=1$, $x+my+z=m$, $x+y+mz=m^2$ :

Solution unique ($m=3$) : resoudre renvoie $[-4,0,1,0,4,0]$, conforme à $(-m-1,1,m+1)$.

Infinité de solutions ($m=1$) : les deux dernières lignes deviennent nulles, aucun pivot n'est trouvé en colonne $2$, et la fonction renvoie None.

Aucune solution : sur $\{\begin{array}{l}x+y=1\\ x+y=2\end{array}$, le même arrêt se produit et la fonction renvoie None.

La limite, à dire honnêtement. La fonction distingue « unique » de « pas unique », mais pas « aucune » de « une infinité » — dans les deux cas elle rend None. Les séparer demande de regarder si une ligne $0=c$ avec $c\ne 0$ subsiste après élimination, ce qui est un ajout de quelques lignes.

Deux choix de calcul. Le pivot partiel (le plus grand coefficient) limite l'amplification des erreurs d'arrondi ; et le test $|\cdot |<\mathtt{tol}$ remplace une comparaison à zéro, illusoire en flottants.