The problem is that the animation variable (t in this case) in animate gets float values, rather than integers.  The simplest way to animate a discrete sequence of plots is to use display rather than animate.  Thus:

> display([seq(PointPlot(L[t]), t=1..6)], insequence=true);

You might try labelledcontourplot from my Maple Advisor Database, <http://www.math.ubc.ca/~israel/advisor>.

After installing the database:

> labelledcontourplot(f, x = 0 .. 4, y = 0 .. 10, filledregions = true, 
       coloring = [yellow, red]);

For a while I thought this was a case of numerical difficulties caused by the very large value of C in your equation, but now I think Maple is right: there really is a singularity here.  I changed epsilon0 to 1E-3 to make things gentler (C is now  79.57747152), and Maple reports the singularity at t = .23214044.  Plotting T(t) for t from 0 to 0.23 shows what looks very much like a vertical asymptote, with T(t) going to -infinity.

You could do it like this:

> L := [a,b,c];
  for v in L do assign(v,Array(1..16)) end do;

Or

> assign(op(L) = seq(Array(1..16),i=1..3));

Puzzle: why does the following not work correctly?

> assign(op(L) = (Array(1..16)$3));

One possible problem I see right away is your line

if  visited[i] != 0 then

"Not equal" in Maple (at least in 1D input)  is <>, not !=

Another problem is your

visited := { op( seq( 0, i = 1..n) ) };

I don't know what you're hoping to accomplish with this: it makes no sense to me.  seq(0, i = 1 .. n) will produce a sequence of n 0's.  op takes either one or two arguments, but produces an error if there are more than two. 

It seems to me that goto can be used within a procedure, but not at top level.  Putting your code in a procedure, I don't get an error.

Your problem can be solved using Laplace transform. 

> PDE:=diff(u(x,t),t$2)+diff(u(x,t),t)-diff(u(x,t),x$2)=sin(Pi*x/l); ICs:=u(x,0)=0,D[2](u)(x,0)=0;
  BCs:=u(0,t)=0,u(l,t)=0;
  with(inttrans):
  Ude:= subs(ICs,laplace(u(x,t),t,s)=U(x),laplace(PDE,t,s));

s^2*U(x)+s*U(x)-diff(U(x),`$`(x,2)) = sin(Pi*x/l)/s

> dsolve({%, U(0) = 0, U(l)=0});

U(x) = l^2*sin(Pi*x/l)/s/((s^2+s)*l^2+Pi^2)

> usol:= invlaplace(rhs(%),s,t);

usol := 1/Pi^2*sin(Pi*x/l)*(l^2-1/(l^2-4*Pi^2)^2*exp(-1/2*t)*((l^2-4*Pi^2)^2*cosh(1/2*t/l^2*(l^2*(l^2-4*Pi^2))^(1/2))*l^2+(l^2*(l^2-4*Pi^2))^(3/2)*sinh(1/2*t/l^2*(l^2*(l^2-4*Pi^2))^(1/2))))

> simplify(usol) assuming l > 2*Pi;

-(4*Pi^2*exp(-1/2*t)*cosh(1/2*t/l*(l^2-4*Pi^2)^(1/2))-4*Pi^2+l^2-exp(-1/2*t)*cosh(1/2*t/l*(l^2-4*Pi^2)^(1/2))*l^2-exp(-1/2*t)*l*(l^2-4*Pi^2)^(1/2)*sinh(1/2*t/l*(l^2-4*Pi^2)^(1/2)))*l^2*sin(Pi*x/l)/(-l^2+4*Pi^2)/Pi^2

Suppose the equation of the plane is a*x + b*y + c*z = d, and your point is [X,Y,Z] with uncertainties dX, dY, dZ (i.e. the x coordinate of the point could be anything from X - dX to X + dX, etc).  Then your point could be in the plane if abs(a*X + b*Y + c*Z - d) <= abs(a)*dX +abs(b)*dY + abs(c)*dZ.

Since jason_cammarata mentioned 126, I suppose this is presumably what he wants, but you might also consider incomplete tictactoe boards
including blanks (which I'll denote by _).  There are 6046 boards in all.

> with(combinat):
  Turns:= X,O,X,O,X,O,X,O,X:
  Boards:= map(op, [seq(permute([Turns[1..j],_$(9-j)]),j=0..9)]):
  nops(Boards);

    6046

Of course not all these boards can arise in an actual game.  Here's a way to check if a board is a win for either player.

> P:= add(B[3*i+1]*B[3*i+2]*B[3*i+3],i=0..2) + add(B[j]*B[3+j]*B[6+j],j=1..3) + B[1]*B[5]*B[9]
    +B[3]*B[5]*B[7];
  IsWin:= proc(L) local LP; LP:= eval(P,B=L); evalb(coeff(LP,X,3)+coeff(LP,O,3)>0) end proc;

Any of the nonwinning boards could arise in a game.

> Wins,NoWins:= selectremove(IsWin, Boards);

A winning board can arise if it comes from a nonwinning board by changing one X or O to a blank.

> OK:= proc(L) ormap(j -> member(subsop(j=_,L),NoWins),[$1..9]) end proc;
  OKWins:= select(OK, Wins):

So here are all the allowed boards.

> OKboards:= [op(OKWins), op(NoWins)]:
  nops(OKboards);

      5478

Manually?  With the mouse.  See the help page ?worksheet,plotinterface,rotateplot

As an animation?  Using the viewpoint option.  See the help page ?plot3d,viewpoint

Maybe something like this?

> kdv:= diff(f(x,t),t) + diff(f(x,t),x,x,x) + 6*f(x,t)*diff(f(x,t),x)=0;
  R:= PDEtools:-TWSolutions(kdv);
  soliton:= subs(R[2], _C1=0, _C2=1, _C3=2, _C4=3, f(x,t));
  plots:-animate(plot, [soliton, x=-4..4], t = -2 .. 2);

Or if you want pdsolve to generate a numerical solution, rather than finding a closed-form soliton:

> S:= pdsolve(kdv,{f(0,t)=f(10,t),D[1](f)(0,t)=D[1](f)(10,t),
    D[1,1](f)(0,t) = D[1,1](f)(10,t),f(x,0) =  cos(x*Pi/5)^2},numeric, spacestep=0.01);
  S:-animate(t=0 .. 5);

What do you mean by "bet that r1 will be negative when r2>0 and r3>0 and that r1 will be positive when r2<0 and r3<0" ?
I take it that r1, r2 and r3 are supposed to be independent random variables, each with probability 1/2 of being positive and 1/2 of being negative, and you win the bet if the outcome is [-,+,+] or [+,-,-] and lose if it is [+,+,+] or [-,-,-].  But what happens when r2 and r3 have different signs?  Do you win then, or is there no bet?

I don't get any error.  No solution either, though.    If you leave out the n[T](0)=0, Maple will "reduce" the problem to a single first-order equation (which may not be very useful to you).  There are probably no "closed-form" solutions.
You could also get series solutions, e.g.:

> dsolve(subs(D=d,{s1,s2,n[T](0)=0}),{n[T](t),n[d](t)},series);

Or:

plots[animate](plot3d, [[r,theta,u(r,t)], r = 0 .. 1, theta = 0 .. 2*Pi, coords = cylindrical], t = 0 .. 1);

What exactly are you trying to do?