The standard method for computing a parametric form of the solution it to use Gaussian elimination for computing a reduced row echelon form of the augmented matrix. Then the unknowns that can be used as parameters are the ones that correspond to columns not containing any leading entry (that is the left most non zero entry in a row or the matrix), and the parametric form can be straightforwardly deduced.
