#My strategy to find a subfield of F_16 with four elements is not too clever. I will create the field F_16 and mod out by my irreducible polynomial, and simply find the element of F_16 whose cubic power is the identity. This will give me a multiplicative group which, along with 0,1 will form the F_4 subgroup.
#So there are some bruteforce themes to my strategy, and I don't want to do all the polynomialmod math by hand. I'll define a quotient ring...
R = PolynomialRing(GF(2),'x')
x = R.gen()
S = R.quotient(x^4 + x + 1, 'a')
a = S.gen()
S

Univariate Quotient Polynomial Ring in a over Finite Field of size 2
with modulus x^4 + x + 1
Univariate Quotient Polynomial Ring in a over Finite Field of size 2 with modulus x^4 + x + 1
