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*}
$$
An unhelpful (?) reduction

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

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.

StackExchange Computer Science Q#131868, answer score: 2