Edit: I've rewritten this whole question because it was unclear. I'm not sure if this is any better. This is a multi-step problem using MATLAB, so I've tried to whittle it down a bit without losing anything.
Let $A$ be a random 6x6 matrix, and $b$ be a random 6x1 matrix. Since $A$ was generated randomly, we would expect it to be nonsingular. The system $Ax = b$ should have a unique solution.
Now let's change $A$ so that it is singular. Replace the third column of $A$ with a linear combination of its first two columns, namely $a_3 = 4*a_1 + 3*a_2$, where $a_1, a_2,$ and $a_3$ are first, second and third columns of $A$.
Now let $y$ be a random 6x1 matrix and let $c = Ay$.
Compute the reduced row echelon form $U$ of $[A$ $c]$.
The free variable determined by the echelon form should be $x_3$. By examining the system corresponding to the matrix $U$, you should be able to determine the solution corresponding to $x_3 = 0$. Let the column vector $w$ be this solution.
Now set the last column of $U$ to zeros. So $U$ should now correspond to the reduced row echelon form of $(A | 0)$. Use $U$ to determine the solution to the homogeneous system when the free variable $x_3 = 1$ and enter your result as vector $z$.
Set $v = w + 3 * z$. The vector $v$ should be a solution to the system $Ax = c$. Why? Explain. . . . What is the value of the free variable $x_3$ for this solution?
Here is the part I'm stuck on:
How could we determine all possible solutions to the system in terms of the vectors $w$ and $z$? Explain.