There really should be a command for this. Explanation follows, but G is your polynomial.
G := Groebner:-FGLM([a^4+a+1, x-a^3], plex(x,a), plex(a,x), characteristic=2);
The problem maps nicely onto an elimination of variables. Given the field equation a^4+a+1=0 and the relation a^3=x, we construct a polynomial system and eliminate a to obtain a polynomial in x only. It just so happens that these equations have a particularly nice form. In general, if you had a big mess you could do the following, and Maple would try to figure out the best way to eliminate variables:
G := Groebner:-Basis([a^4+a+1, x-a^3], plex(a,x), characteristic=2);
In this case though you have a triangular system of equations in which each variable appears as the largest term of one equation. This is a lex Groebner basis with x > a. The command at the top applies a linear algebra method for converting Groebner bases, literally searching for linear dependencies among the powers of x=a^3. It finds one at x^4.
It uses sparse linear algebra. It should be reasonably fast.