Algorithm that solves a system of linear equations over finite fields when a parameter is needed

I was reading Kipnis' and Shamir's paper on Cryptanalysis of the HFE Public Key Cryptosystem by Relinearisation and I wanted to implement the example at the end in Octave without using any additional packages (e.g symbolic ...). I would like to create an algorithm that solves systems of linear equations over finite fields (q = 7 in this case) where you have more variables than equations (in this case I am going to need one parameter).

I'm pretty new to this topic so I tried setting $y_{12} = z$ (example from the paper) and subtracted the $y_{12}$ vector from the solution for $z = ${$1, 2$} and solved the system of equations for the two new solutions to try to find the slope of the parametric solutions:
$ y_{11} = 2 + 5z \\ y_{12} = z \\ y_{13} = 3 + 2z \\ ... $

This didn't seem to work for reasons which I think I know and now I'm clueless on what to do. I'd appreciate your help.
To save you a bit of time copying the equations onto Octave, I'll leave you with this:

eq = [3 5 5 2 6 4 5; 6 1 4 4 5 1 6; 5 2 6 2 3 2 5; 2 0 1 6 5 5 0; 4 6 2 5 1 4 0];