This is a system of linear equations, so you can write it as $A(a_0,a_1,a_2) = b$
var('x')
F.<x> = GF(2^8, modulus=x^8 + x^5 + x^3 + x + 1)
R.<a0,a1,a2> = PolynomialRing(F)
eqns = [a0+a1*x+a2*x^2 - (x^7+x^6+x^5+x^2+x), a0+a1*x^2+a2*x^4 - (x^7+x^6+x^5+x^2), a0+a1*x^3+a2*x^6 - (x^7+x^3+x^2+1)]
A = matrix(F, [[eqn.coefficient(b) for b in R.gens()] for eqn in eqns])
b = vector(F, [-eqn.constant_coefficient() for eqn in eqns])
and then you can solve it:
sage: A.solve_right(b)
(x^7 + x^2 + 1, x^7 + x^2 + 1, x^7 + x^2 + 1)
You can also do the following, because the ideal that these equations define is zero-dimensional:
sage: var('x')
x
sage: F.<x> = GF(2^8, modulus=x^8 + x^5 + x^3 + x + 1)
sage: R.<a0,a1,a2> = PolynomialRing(F)
sage: I = R.ideal([a0+a1*x+a2*x^2 - (x^7+x^6+x^5+x^2+x), a0+a1*x^2+a2*x^4 - (x^7+x^6+x^5+x^2), a0+a1*x^3+a2*x^6 - (x^7+x^3+x^2+1)])
sage: I.dimension()
0
sage: I.variety()
[{a1: x^7 + x^2 + 1, a2: x^7 + x^2 + 1, a0: x^7 + x^2 + 1}]