A trick that can be used here is to use the fact that the polynomial
$$L(x):=x^2+(\alpha^3+\alpha)x$$
has the property
$$L(x+y)=L(x)+L(y)$$
for all $x,y\in F$. This is because in characteristic two fields we have the rule
$(x+y)^2=x^2+y^2$.
So $L(x)$ is a linear transformation (over the prime field $\Bbb{Z}_2$). We can use tools from linear algebra. Let's check what $L$ does to the basis elements:
$$
\begin{aligned}
L(1)&=\alpha^3+\alpha+1,\\
L(\alpha)&=\alpha^4=\alpha+1,\\
L(\alpha^2)&=\alpha^5+\alpha^4+\alpha^3=\alpha^3+\alpha^2+1,\\
L(\alpha^3)&=\alpha^4=\alpha+1.
\end{aligned}.
$$
We pretty much see right away that
$$L(\alpha+\alpha^2)=L(\alpha^2)+L(\alpha)=\alpha^3+\alpha^2+\alpha$$
(we could also do this with linear algebra and coordinates as was implied in Lonza Leggiera's answer, +1) giving us one root $x_1=\alpha^2+\alpha$.
Obviously $L(\alpha^3+\alpha)=(\alpha^3+\alpha)^2+(\alpha^3+\alpha)^2=0$. Because $L$ is a quadratic, it can have at most two zeros, so the kernel of $L$ is 1-dimensional, spanned by $\alpha^3+\alpha$.
This means that the other solution is
$$
x_2=x_1+(\alpha^3+\alpha)=\alpha^3+\alpha^2.
$$