Thank you very much for the help. This works like a charm.

Thank you very much for the procedure! I am using LinearAlgebra[Generic] since I am working over finite fields. Somehow, the procedure does not seem to work (although in the standard LinearAlgebra package it works) # Example in LinearAlgebra that works with(LinearAlgebra): nestedseq:= proc(numtimes::nonnegint, middle) if numtimes = 0 then middle else seq(eval(nestedseq(numtimes-1, middle)), cat('j',numtimes)=[1,2]) end if end proc; nestedseq(4, Determinant(Matrix([[j1,j2],[j3,j4]]))); #Example in LinearAlgebra[Generic]: modify nestedseq: nestedseq:= proc(numtimes::nonnegint, middle,felems) if numtimes = 0 then middle else seq(eval(nestedseq(numtimes-1, middle,felems)), cat('j',numtimes)=felems) end if end proc; # Then we do the preliminaries to build GF4 >restart; with(LinearAlgebra[Generic]): p := 2; m4 := z^2+z+1; felemGF4:=[0,1,z,z+1]; GF4[`0`] := 0: GF4[`1`] := 1: GF4[`=`] := `=`: GF4[`+`] := ()->modp(`+`(args),p): GF4[`-`] := (a,b)->`if`(nargs=1,modp(-a,p),modp(a-b,p)): GF4[`*`] := (a,b)->modp(Rem(a*b,m4,z),p): GF4[`/`] := proc(a,b) local i; if b=0 then error "division by zero"; end if; Gcdex(b,m4,z,'i') mod p; Rem(a*i,m4,z) mod p; end proc: # Say we calculate some determinants >nestedseq(4, Determinant[GF4](Matrix([[j1,j2],[j3,j4]])),[0,1,z,z+1]); gives the following error: Error, (in mod/Expand) polynomial over a finite field expected