For your first computation, Maple does give me 1, as shown by the following code:

restart:
p:= 3:
firred:= T^2+2*T+2:
d:= degree(firred):
q:= p^d:
alias(tau = RootOf(firred)):
Expand((tau+1)*(2*tau+2)) mod 3;

The Expand could also be Normal in this example. 

Your error is caused by using 2D Input (not your fault). In 2D Input, (tau+1)(2tau+2) is not interpretted as a product of two binomials. Unfortunately, it's interpretted as a "function" (tau+1) being applied to an argument (2tau+2). The degenerate "function" is a constant function, so the result will be tau+1 regardless of what's in the second pair of parentheses. There are two ways to correct this:

1. Using 2D Input, put a space between the two pairs of parentheses: (tau+1) (2tau+2). In all cases where there's an ambiguity between implied multiplication and function application, a space is needed to force it to be implied multiplication (no matter how stupid or degenerate the "function" may be). For example 2(x+3) needs to be 2 (x+3).
2. In any form of input, you can use * for multiplication. You just need to get out of the habit of using juxtaposition as implied multiplication.

Regarding your second computation, I have done it meticulously with pen and paper, and I get what Maple gets: Rem(x^q, f, x) mod 3 is x, not (tau+1)*x. You'll need to show me the computation that leads you to believe your answer. Indeed, my pen-and-paper computation confirms all 9 entries in this table:

f:= x^2 + tau*x:
interface(rtablesize= 2+q):
<
   <` k ` | x^` k`>,
   <`===` | `=====`>,
   < 
      <seq(1..q)> | 
      Vector(q, k-> collect(Rem(x^k, f, x) mod p, x))
   >
>;

Using anything may "work"---if by "work" you simply mean "returns true"---but that's not achieving the effect that you want, so I wouldn't say that it "works correctly". The coefficients of g aren't simply tau; they are themselves polynomials in tau with integer coefficients. So, what you want is

type(g,  polynom(polynom(integer, tau), x));

I'm not sure exactly what you mean, but I think that what you want is given by the command simplex[ratio]. Certainly, this command returns a list of multipliers. See ?simplex and ?simplex,ratio.

The command select can be used to link your two pieces of code. The general form is select(P, L), where P is a boolean procedure (i.e., one that returns true or false) acting on the elements of a container L (such as a list or set). 

I wouldn't use your code above, because both pieces are extremely inefficient (although they do perform the computation correctly, which is the most important thing, and they're both fairly clean stylistically), but if I was going to use it with select, I'd do

select(p-> findOrderOf(2,p) = p-1 and findOrderOf(3,p) = p-1, primes),

which can be improved with andmap to 

select(p-> andmap(g-> findOrderOf(g,p) = p-1, [2,3]), primes).

Your order-finding procedure is inefficient for a mathematical reason: It doesn't make use of the fact that the order must be a divisor of the totient^[*1] of the modulus, which means that there are only a few exponents that actually need to be checked.

Your prime-listing code is inefficient due to an idiosyncracy of Maple rather than a mathematical reason: Lists (and sequences and sets) are effectively read-only; they cannot be added to (or even changed in any way). Any attempt to do it actually creates a fresh list, copying the parts of the original that are needed, no matter how many elements need to be copied. That's done in the background, transparently to the user, who often doesn't realize that something rather inefficient is occuring. Some other Maple containers such as Vectors and tables do not have this read-only behavior.^[*2]

So here's some efficient code: If you're simply interested in answering the given homework problem, you could do

select(
   p-> isprime(p) and NumberTheory:-MultiplicativeOrder~({2,3}, p) = {p-1}, 
   [seq(2001..2999, 2)] #only need to check odds
);  

If you'd rather not rely on Maple's prewritten MultiplicativeOrder, you could do

IsGenerator:= (x::integer, p::prime)->
   not ormap((d,g)-> g &^ d mod p in {0,1}, (p-1) /~ op~(1, ifactors(p-1)[2]), irem(x,p))
:
select(p-> isprime(p) and andmap(g-> IsGenerator(g,p), [2,3]), [seq(2001..2999, 2)]);

My procedure IsGenerator makes use of a special efficiency that applies only when the modulus p is prime: that the totient in that case is simply p-1. There's also an efficiency (regardless of whether the modulus is prime) that results from if the order is not the totient, then we know that x is not a generator without needing to actually compute its order.

The special operator &^ (as in g &^ d mod p) computes g^d mod p without first computing g^d. Thus, it can be used efficiently even when the exponent is extremely large. (The operator itself is inert: g &^ d by itself, without mod, does nothing.)

^[*1]As you probably know, the totient of a positive integer n is the count of positive integers x <= n such that gcd(x,n) = 1, i.e.,  x and n are relatively prime. Thus, the totient of n is the order of the multiplicative group mod n. The totient is also called "Euler's phi function".

^[*2]In Maple documentation, structures that are not read-only (in the way that I described above) are called mutable, and the others are called immutable (or perhaps non-mutable). I think that the term read-only is more widely known and the concept more widely understood.

There are two problems with the worksheet that you posted, one typographical and one logical. For the typographical problem, look closely at the line that begins consd1:=. Notice how that line is printed in a different way than the line immediately above it? I think that that means that that line is considered textual commentary, and so it isn't being executed. I suggest that you simply retype that line. You can ensure that it gets placed in an execution group by using Ctrl-J (while your cursor is on that line) and retyping it at the prompt that appears.

The logical problem is that the variables x[2,1], x[3,1], x[4,1], and x[5,1] do not have any constraints. Hence they can take any negative values, making the problem unbounded towards negative infinity.

Your line

cons[k] := add(x[i, k], i = N)+add(x[k, j], N) = 2

should be

cons[k] := add(x[i, k], i = N)+add(x[k, j], i= N) = 2​​​​​​​

As a workaround, the original 33x32 system can be easily and quickly solved by the direct use of the LinearAlgebra package:

macro(LA= LinearAlgebra):
(A,B):= LA:-GenerateMatrix(liseq1, lisvar1):
(P,L,U):= LA:-LUDecomposition(A):
Sol:= lisvar1 =~ seq(LA:-BackwardSubstitute(U, LA:-ForwardSubstitute(L, P^+.B))):
andmap(evalb, eval(liseq1, Sol)); #verify correctness

A small adjustment to Kitonum's second method, the one that uses subsindets, will make it safer (more robust) in that the conversion of  `*` to `.` will only occur when intended:

subsindets['flat'](A.B, 'specop'(rtable, `*`), `.`@op);

If your Vectors have purely numeric entries, this won't make any difference; but if they themsleves contain products, it might.

Maple's assignment operator is two ordinary ASCII characters: a colon (:) followed by an equals (=). The character that you have after X is some exotic non-ASCII character that looks like :=, but isn't.

If you make this small change, then your example will work as expected.

Subscripting is not quite the same thing as indexing, and the underscore _ is not an indexing operator in Maple (nor is it a subscripting operator, though double underscore __ is). 

There are many ways to do what you want:

Most-straightforward way:

seq([X[k], Y[k]], k= 1..nops(Y));

Clever, minimal-code, way:

`[]`~(X,Y)[]:

Perhaps the most efficient (although I haven't checked):

convert(<<X> | <Y>>, list, 'nested')[];

Most versatile: Just make it into a two-column Matrix:

XY:= <<X> | <Y>>:

This last form is accepted by most Statistics commands and plot; indeed, they are more efficient when the input data is presented thus.

As you can see, of the four methods, only one uses an index. I tend to prefer the clean syntax of these index-free methods over those that require introducing a bound variable for the index; however, these latter methods are sometimes a little more efficient.

One way is to use Rem with mod. Here's an example over GF(3,2):

restart:
alias(alpha= RootOf(Randprime(2,z) mod 3)):
f:= x*(x-1):
Rem(alpha^5, f, x) mod 3;

                          2*alpha

In addition, you could redefine your ring operations `+`, `*`, `-` to perform this Rem operation; or you could give them other names, which would still be simpler than frequently typing Rem(...., f, x) mod 3. If the ring has only a few hundred elements or fewer, remember tables for the operations are very effective for high-speed computations. 

My preference is 

n:= 10:
cat(seq("TH"[x], x= rtable(1..n, random(1..2))));

I realize that this probably looks weird to most readers. I like it because I have a suspicion (not confirmed yet) that rtable with a random or frandom initializer is very efficient, more efficient than a sequence of calls to a procedure produced by rand.

The following code is equivalent to your code except that the recursive procedure GoDownTree returns a sequence. Therefore, each level of recursion doubles the length of that sequence. You are correct that print is never a good way for any procedure to return its value, let alone a recursive procedure.

istZweirpotenz:= (n::posint)-> evalb(n = Scale2(1, ilog2(n))):

SubProdTree:= proc(n::satisfies(istZweirpotenz), x::name, u::And(list, satisfies(u-> nops(u)=n)))
local k:= ilog2(n), j, i, M::table(polynom(integer, x)), leftList, rightList, kList;
   for j to n do M[0, j-1]:= x-u[j] od;
   for kList from 0 to k do
       leftList[0, kList]:= M[0, 2*kList];
       rightList[0, kList]:= M[0, 2*kList+1]
   od;
   for i to k do
      for j from 0 to 2^(k-i)-1 do M[i,j]:= M[i-1, 2*j]*M[i-1, 2*j+1] od;
      for kList from 0 to 2^(k-i-1)-1 do
         leftList[i, kList]:= [leftList[i-1, 2*kList], rightList[i-1, 2*kList], M[i, 2*kList]];
         rightList[i, kList]:= [leftList[i-1, 2*kList+1], rightList[i-1, 2*kList+1], M[i, 2*kList+1]]
      od
   od;
   ListTools:-Flatten([leftList[k-1, 0], rightList[k-1, 0], M[k,0]])
end proc:         

GoDownTree:= (f::polynom, n::satisfies(istZweirpotenz), x::name, l::list, p::prime)->
   if n = 1 then f 
   else 
      thisproc(evala(Rem(f, l[n-1], x) mod p), n/2, x, l[..n-1], p),
      thisproc(evala(Rem(f, l[2*n-2], x) mod p), n/2, x, l[n..2*n-2], p)
   fi
:
FastEval:= (f::polynom, n::satisfies(istZweirpotenz), x::name, p::prime, u::list)->
   GoDownTree(f, n, x, SubProdTree(n, x, u), p)
:

FastEval(x^2+x+1, 8, x, 11, [$1..8]);

      3, 7, 2, 10, 9, 10, 2, 7

I note that your code did not make any use of the field GF(11,2). So, my code doesn't either, as I said my code is equivalent to yours. So, the computations are done over GF(11,1), but I don't know if that was your intention. The evala command does nothing in this example. As far as I know, there is never a reason to use both evala and mod, but I am not sure about that.

Pay close attention to the ways that I changed your code, especially:

1. ilog2(n) is equivalent to floor(log[2](n)) when n is a positive integer, but it is many times faster.
2. Procedures usually should give an error when they are given invalid input, not return false.
3. Recursive procedures should almost always return values.
4. Unassigned global names (x in this case) should be passed in rather than having the procedure just assume that they're unassigned.
5. Do not use the with command to load packages for use in procedures.

The issue is the lack of a space character after &t, as pointed out by Torre. Whether you use a space before &t is unimportant.

The code layout on those online help pages looks like some automatically generated crap, and it's easy to see how they'd lead you to make such an error. I recommend that you always compare those pages with the in-program help, in this case ?DifferentialGeometry,Tensor,KillingVectors. For additional clarity, press F5 so that the code of the help page examples will be displayed in a monospaced font (if it isn't already thus displayed).

I'm not aware of any situation in Maple where an uninterupted sequence of alphabetic letters is interpretted as multiple symbols. So, the way that you've coded it, Maple isn't even aware that you're trying to use &t as a neutral operator; rather, it thinks that you're trying to use &tdx as a neutral operator, which is invalid because it's not followed by an operand.

@mmcdara Yes, what you ask for is essentially possible if you're willing to accept an equivalent of "indeterminate" instead of undefined.

solve({a*b = a*c}, b, parametric);

      piecewise(a = 0, [{b = b}], a <> 0, [{b = c}])

It's very rare (if at all) that much useful can be done in Maple by assuming a one-point inequality. And solve is quite resistant to all assumptions, although there are sometimes ways around that.