Complexity of a decision problem: system of linear equations over finite field with restricted solutions

I have a system of linear equations over a finite field $$\mathbb F_p \cong \mathbb Z_p$$, and I'm interested in the decision problem of whether there exists a solution where all of the variables $$x_i$$ are in the set $$\{0, 1\} \subset \mathbb F$$. In particular, I'm trying to determine whether this problem is $$\mathcal{NP}$$-hard.

Example

One system of equations over $$\mathbb F_3$$ is: \begin{alignat*}{2} &x_1\begin{bmatrix}1 \\ 1 \\ 0\end{bmatrix} + x_2\begin{bmatrix}0 \\ 0 \\ 0\end{bmatrix} + x_3\begin{bmatrix}0 \\ 2 \\ 0\end{bmatrix} + x_4\begin{bmatrix}1 \\ 0 \\ 1\end{bmatrix} \\ + \, &x_5\begin{bmatrix}2 \\ 1 \\ 0\end{bmatrix} + x_6\begin{bmatrix}0 \\ 1 \\ 1\end{bmatrix} + x_7\begin{bmatrix}1 \\ 0 \\ 0\end{bmatrix} + x_8\begin{bmatrix}1 \\ 2 \\ 0\end{bmatrix} = \begin{bmatrix}0 \\ 0 \\ 2\end{bmatrix}. \end{alignat*} This system of equations is satisfiable with entries in $$\{0,1\}^8 \subset \mathbb F^8$$, namely \begin{align*} (x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8) &= (0,0,0,1,0,1,1,1) \hspace{1em}\text{or}\\ (x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8) &= (0,1,0,1,0,1,1,1). \end{align*}

One suggestion that was given to me was turning this into a system of quadratic equations in the following way: define auxiliary functions coordinatewise \begin{align*} w_1 &= x_1 + x_4 + 2x_5 + x_7 + x_8 \\ w_2 &= x_1 + 2x_3 + x_5 + x_6 + 2x_8 \\ w_3 &= x_4 + x_6 - 2, \end{align*} and use these to solve the system of quadratic and linear equations $$w_1 = w_2 = w_3 = x_1^2 - x_1 = \cdots = x_8^2 - x_8 = 0.$$

However, the MQ-problem (Multivariate Quadratic equations over a finite field) is $$\mathcal{NP}$$-hard, so this reduction doesn't help. However, this set-up is a quite special case, so I'm holding out some hope that the original problem might still be in $$\mathcal{P}$$.

Is there a polynomial-time algorithm for determining the existence of a solution of linear equations over a finite field with restricted variables? Or is it known if this problem is $$\mathcal{NP}$$-hard like the MQ-problem?

You can reduce 1-IN-3 SAT to your problem (an instance is a 3CNF, and we want to find a satisfying assignment having exactly one satisfied literal per clause), assuming $$p \geq 3$$.
A clause $$x \lor y \lor z$$ is encoded as the constraint $$x+y+z=1$$.
A clause $$\bar x \lor y \lor z$$ is encoded as the constraint $$1-x+y+z = 1$$; and so on.
When $$p = 2$$, your problem becomes easy.