I am working through some elementary linear algebra exercises and I want to obtain the solution set to the following system of linear equations using MatLab’s solve() function, from the symbolic toolbox add-on:

2x      + z + w = 5
     y      - w = -1
3x      - z - w = 0
4x + y + 2z + w = 9

This system, when converted to reduced row-echelon form, retains a free variable, meaning there are infinitely many solutions, but when I use solve(), it simply returns a particular solution to the equation:

syms x y z w h i j k
eq1 = 2*x + z + w == 5;
eq2 = y - w == -1;
eq3 = 3*x - z - w == 0;
eq4 = 4*x + y + 2*z + w == 9;
solve([eq1,eq2,eq3,eq4])

returns

struct with fields:
    w: 3
    x: 1
    y: 2
    z: 0

This is a solution, but not the only solution. I would like to be able to (a) obtain the full solution set using solve(), and (b) not have to check whether the system has infinitely-many solutions outside of using solve(). (Currently, I am also setting up an augmented coefficient matrix, A, and using rref(A) to see whether there are any free variables). Are these goals possible or should I perhaps be using another function?