Contradiction when solving a linear system with Gauss-Jordan elimination Consider a linear system with an unknown constant $a$:
$$
\left\{
\begin{aligned}
x+(a-1)y+az&=1 \\
ax+ay+az&=1 \\ 
a^2x+y+z&=a
\end{aligned}
\right. 
$$
This gives us an augmented matrix:
$$A= \left[
\begin{array}{ccc|c}
  1&a-1&a&1\\
  a&a&a&1\\
  a^2&1&1&a
\end{array}
\right] $$
From the augmented matrix, it is clear to me that when $a = 1$, we will have one parameter, since there are three unknowns but only two equations.
Let $a = 1$:
$$ \left[
\begin{array}{ccc|c}
  1&0&1&1\\
  1&1&1&1\\
  1&1&1&1
\end{array}
\right] \xrightarrow{R_3 - R_2}
\left[
\begin{array}{ccc|c}
  1&0&1&1\\
  1&1&1&1\\
  0&0&0&0
\end{array}
\right] \xrightarrow{R_2 - R_1}
\left[
\begin{array}{ccc|c}
  1&0&1&1\\
  0&1&0&0\\
  0&0&0&0
\end{array}
\right]
 $$
So we get the solution $z = s, y = 0, x = 1 - s \implies \text{Infinite number of solutions}$.
However, if instead I don't consider $a = 1$ but only $a \ne 0$, then I can proceed with Gauss-Jordan elimination to eventually get the following RREF after these series of elementary row operations:
$$
\frac1aR_2\longrightarrow
R_1 \leftrightarrow R_2 \longrightarrow 
R_2-R_1 \longrightarrow 
R_1-\frac{1}{a^2}R_3 \longrightarrow 
\frac{1}{a^2}R_3 \longrightarrow 
R_3+R_1 \longrightarrow 
a^2R_1 \longrightarrow 
\frac{1}{a^2-1}R_1 \longrightarrow 
R_2-aR_1 \longrightarrow 
R_2+2R_1 \longrightarrow 
R_3\leftrightarrow 
R_1 \longrightarrow 
R_2\leftrightarrow 
R_1 \longrightarrow 
R_3-R_2 \longrightarrow 
R_2-R_1
$$
$$A_\text{RREF}=
\left[
\begin{array}{ccc|c}
  1&0&0&\frac1a\\
  0&1&0&\frac{-(a-1)}{a}\\
  0&0&1&\frac{a-1}{a}
\end{array}
\right]
$$
I had assumed that there was an error in my Gauss-Jordan elimination because this RREF implies that when $a = 1$, we have one and only one unique solution of $x = 1, y = 0, z = 0$. However, I plugged this matrix into MatLab and got the same RREF.
How can it be that when $a = 1$, we get two different solutions from the two row-equivalent matrices? 

If you want to rref and avoid any division at all you'll get
$$\begin{pmatrix}
a(a-2)(a^2-1)& 0 & 0 & (a-2)(a^2-1)\\
0 & -a(a-2)(a+1) & 0 & (a-2)(a^2-1)\\
0 & 0 & a(a^2-1) & (a+1)(a-1)^2
\end{pmatrix}.$$
From here it's eazy to consider all different cases.
