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?

A correlation coefficient does not determine the joint distribution of the coins.  I'll assume that this distribution is symmetric with respect to permutations of the coins.  Then there are three degrees of freedom: we may take P(HHH) = a, P(HHT) = P(HTH) = P(THH) = b, P(HTT) = P(THT) = P(TTH) = c, P(TTT) = d where a,b,c,d >= 0,
a + 3*b + 3*c + d = 1.  If X_i are the Bernoulli random variables corresponding to heads for coins 1,2,3, then the correlation coefficient of any two of these is
 

corr := (a+b-(a+2*b+c)^2)/(a+2*b+c-(a+2*b+c)^2)

You could specify the correlation coefficient and the expected value (i.e. the probability of heads for each coin), and still have one degree of freedom.
Thus with correlation coefficient p and expected value m:
 

> eqs:= {a+3*b + 3*c + d = 1, a + 2*b + c = m, corr = p};
  solve(eqs,{a,b,c});

{a = -3*m-d+1+3*m^2+3*p*m-3*p*m^2, b = 3*m+d-1-2*m^2-2*p*m+2*p*m^2, c = -2*m-d+1+m^2+p*m-p*m^2}

The probability of three heads or three tails is then a + d = -3*m+1+3*m^2+3*p*m-3*p*m^2.