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? 
 A: The contradiction here, if any, is between your expectation of what MatLab does and what it actually does. You are correct that there are multiple cases to be considered, but MatLab and other software plunges ahead without taking potential divisions by zero into account when the matrix elements involve unbound variables. From the MatLab documentation for rref: 

If the elements of a matrix contain free symbolic variables, rref regards the matrix as nonzero.

In fact there are more cases to consider than the two you mentioned. If $a=0$, then the second RREF in your question involves division by zero, so you need to examine that case separately, too.
A: In the case $a=1$, the rows 2 and 3 of the $3\times 3$ matrix are the same and so the matrix has rank 2 and there are infinitely many solutions (the solution set is a line).
If $a\ne 1$, the $3\times 3$ matrix has rank 3 and there is unique solution. 
A: 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.
