We would like to find a matrix $\mathrm X \in \mathbb F_2^{4 \times 4}$ such that $\mathrm X 1_4 = 1_4$ and $1_4^{\top} \mathrm X = 1_4^{\top}$. Vectorizing, we obtain a system of $8$ linear equations in $16$ unknowns
$$\begin{bmatrix} 1_4^\top \otimes \mathrm I_4\\ \mathrm I_4 \otimes 1_4^\top\end{bmatrix} \mbox{vec} (\mathrm X) = \begin{bmatrix} 1_4\\ 1_4\end{bmatrix}$$
Using Gaussian elimination (over $\mathbb F_2$), we eventually obtain the augmented matrix in RREF
$$\left[\begin{array}{cccccccccccccccc|c} \boxed{1} & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 0\\ 0 & \boxed{1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1\\ 0 & 0 & \boxed{1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 1\\ 0 & 0 & 0 & \boxed{1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 0 & \boxed{1} & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boxed{1} & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boxed{1} & 1 & 1 & 1 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$
Visual inspection of the augmented matrix in RREF allows us to conclude that the linear system is consistent, that none of the columns corresponding to non-pivot variables are all-zero and that
$$\mbox{rank} \begin{bmatrix} 1_4^\top \otimes \mathrm I_4\\ \mathrm I_4 \otimes 1_4^\top\end{bmatrix} = 7$$
Thus, we have $16 - 7 = 9$ degrees of freedom and $\color{blue}{2^9 = 512}$ solutions.
Addendum
Adding the 1st, 2nd, 5th and 6th columns of the augmented matrix in RREF produces the zero vector. Adding the 1st, 3rd, 5th and 7th columns of the augmented matrix in RREF also produces the zero vector. Proceeding in this manner, we can easily find a basis for the $9$-dimensional null space.
Unvectorizing each basis vector, we obtain the $9$ basis matrices
$$\mathrm N_1 := \begin{bmatrix} 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_2 := \begin{bmatrix} 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_3 := \begin{bmatrix} 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\end{bmatrix}$$
$$\mathrm N_4 := \begin{bmatrix} 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_5 := \begin{bmatrix} 1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_6 := \begin{bmatrix} 1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\end{bmatrix}$$
$$\mathrm N_7 := \begin{bmatrix} 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_8 := \begin{bmatrix} 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\end{bmatrix} \qquad \qquad \mathrm N_9 := \begin{bmatrix} 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\end{bmatrix}$$
Hence, the set of all $2^9 = 512$ solutions of $\mathrm X 1_4 = 1_4$ and $1_4^{\top} \mathrm X = 1_4^{\top}$ can be generated as follows
$$\left\{ \mathrm X_p + \sum_{k=1}^9 z_k \mathrm N_k \,\, \bigg| \,\, z_k \in \{0,1\} \right\}$$
where $\mathrm X_p$ is a particular solution, e.g.,
$$\mathrm X_p := \begin{bmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 1 & 1\end{bmatrix}$$