To construct GF(16) as a two-step extension, which is what it appears you are interested in, would go as follows.

```
sage: k=GF(2)
sage: kx.<x>=k[]
sage: kt.<t>=k[]
sage: kx.<x>=k.extension(t^2+t+1)
sage: x^2+x+1
0
sage: kxu.<u>=kx[]
sage: kxy.<y>=kx.extension(u^2+x*u+1)
sage: y^2+x*y+1
0
sage: kxy
Univariate Quotient Polynomial Ring in y over Finite Field in x of size 2^2 with modulus y^2 + x*y + 1
```

This has a slight problem: the field kxy constructed above is a relatively inefficient version of GF(16): sage just knows it's a quotient of a polynomial ring, whereas for GF(16) sage has a specialized implementation that should be quite fast.

As suggested by dan_fuela one you can just construct GF(16) any way you want and then find the elements x,y with the desired properties in it. You don't have to exhaustively search the whole field for the right minimal polynomials to occur though: you can just ask sage to solve the relevant polynomial equation. In GF(16) it doesn't matter much because it's a tiny field, but enumerating larger fields quickly becomes undoable.

```
sage: kxy=GF(16)
sage: kxyT.<T>=kxy[]
sage: x=(T^2+T+1).roots()[0][0] # this just selects the first root from the full list of roots with multiplicities
sage: y=(T^2+x*T+1).roots()[0][0]
sage: x^2+x+1
0
sage: y^2+x*y+1
0
```