You can use Gaussian elimination. Write it is a matrix equation - but it is a 3x4 matrix, 3 rows and 4 columns.
Subtract two times the first row from the second and add the first to the third. The gives you just a 1 in the first column.
You will then find that the (2,2) entry becomes 5 and the (3,2) becomes 1. Divide row 2 by 5 to get a 1 in position (2,2); now subtract row 2 from row 1 and row 3. So far we have a 1 at the top of the first column, then 0s, and a 1 in the middle of the second column and 0s above and below.
Do this again for the third column.
This gives the identity matrix on the left, followed by a column of three numbers, say $a$, $b$ and $c$; suppose the right hand side becomes $A$, $B$ and $C$. We then have the general solution in the form
$$(x_1, x_2, x_3, x_4) = (A, B, C, 0) + \lambda (a, b, c, -1).$$
Note that we could have done this differently, say leaving the 2nd column and 'clearly' the other three; this would give a different, but equivalent, form of the answer. This is because your solution is a line of solutions. Doing it another way, you will always get the ratio $a : b : c : -1$ the same, but you may get a different factor.
By considering the initial work in crash's answer, the solution is
$$\left(\begin{matrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{matrix}\right) = \left(\begin{matrix} 1 \\ 0 \\ 0 \\ 0 \end{matrix}\right) + \lambda\left(\begin{matrix} -1 \\ 2 \\ -1 \\ 0 \end{matrix}\right) + \mu\left(\begin{matrix} -1 \\ 2 \\ 0 \\ -1 \end{matrix}\right),$$
where varying $\lambda$ and $\mu$ give you your two dimensional solution space.
As crash said, this is one of those things that it is worth showing at least once in your life, but then you just use it.
Note that in general, given 4 variables and 3 equations, you will have a 1-dimensional solution space. In this case it is 2D because any 3x3 submatrix has determinant zero, but there exist 2x2 submatrices with non-zero determinant.
Hopefully this helps! If it does, then please remember to upvote and/or accept this answer! :)