Utilizamos el teorema Rouch-Frobenius, de modo que calculamos el determinante 
de la matriz del sistema y para ello la introducimos y calculamos el determinante

matrix(
 [1,1,a], 
 [1,a,1], 
 [a,1,1]
);

determinant(%);

Vamos a ver cuando este determinante es cero; y para eso igualamos a cero y resolvemos
(recordemos que la variable es a)

solve([%o2], [a]);

El determinante no es cero si a es distitnto de -2 y distinto de 1. Por tanto, 
en ese caso el sistema es compatible determinado ya que el rango de la matriz es 3 y
el de la ampliada también. Podemos entonces resolver el sistema para estos casos.

linsolve([x+y+a*z=a^2, x+a*y+z=a, a*x+y+z=1], [x,y,z]);

Vamos a estudia el caso a=-2, entonces el sistema es x+y-2*z=4; x-2*y+z=-2; -2x+y+z=1. 
Observemos que, como el menor formado por las dos primeras filas y las dos primeras 
columnas es distinto de cero (veámoslo)

matrix([1,1],[1,-2]);

determinant(%);

POdemos sustituir la ultima columna de la matriz del sistema por la columna de 
términos independientes y estudiar su determinante

matrix(
 [1,1,4], 
 [1,-2,-2], 
 [-2,1,1]
);

determinant(%);

Como este determinante es distitnto de cero, el rango de la matriz ampliada 
es 3 y por tanto, si a=-2, el sistema es incompatible.

Vamos a estudiar ahora el caso a=1. El sistema tiene todas sus ecuaciones igales a 
x+y+z=1, con lo que es claramente indeterminado. Lo resolvemos

linsolve([x+y+z=1], [x,y,z]);

Que, como era de esparar depende de dos parámetro %r1 y %r2