> de:= 4*(D@@2)(x)(t)+3*D(x)(t)+k*x(t) = 0;
  characteristicpoly:= normal(eval(lhs(de), x = (t -> exp(r*t)))/exp(r*t));
  underdamped:= solve(discrim(characteristicpoly,r) < 0, {k});

This is the condition for the equation to be underdamped.  OK, take a suitable k, say 2.

> sol := dsolve({eval(de, k=2),x(0)=0,D(x)(0)=1});

sol := x(t) = 8/23*23^(1/2)*exp(-3/8*t)*sin(1/8*23^(1/2)*t)

  Now surely you can see that this will be 0 infinitely many times. But if you want Maple to find the zeros:

> solve(rhs(sol) = 0, t, AllSolutions);

Note that a parameter such as _Z1~ stands for "an arbitrary integer".

2D math input strikes again!   You have ρsl with no spaces, which Maple treats as a single variable.  This is one reason why I strongly recommend using 1D maple input. 

Also there is confusion between

   p S U
   -----
     2
 ---------
      m
--------------
   p s l
   -----
     2

and

    p S U
    -----
      2
-------------
      m
   --------
     p s l
     -----
       2

For a Matrix with a full basis of eigenvectors (which this one will have with high probability) you can proceed as follows:

> E,V := Eigenvectors(M):
  V . DiagonalMatrix(map(e -> exp(-I*t*e),E)) . V^(-1);

This should be very fast.  It may be less accurate in some cases than using MatrixExponential, though.

A zero-sum game corresponds to a linear programming problem. 

There are four pure strategies for Player 1, which I'll denote [a,b] where a and b are either 0 = Play or 1 = Pass (a for heads and b for tails). 
The pure strategies for Player 2 are 0 (play) or 1 (pass). 
When Player 1 plays [a,b] and Player 2 plays c, if the coin is heads the payoff for Player 1 is 

Pheads:= - a + (1-a)*c + 2*(1-a)*(1-c);

If it is tails the payoff for Player 1 is

Ptails:= -b  + (1-b)*c - 2*(1-b)*(1-c); 

Since it is a fair coin, the expected payoff in this case is 1/2*Pheads + 1/2*Ptails.

Now the linear programming problem for Player 1 is to find a mixed strategy to maximize the minimum payoff over all of Player 2's pure strategies.
Thus:

constraints:= {seq(v <= add(add(x[a,b]*(Pheads + Ptails)/2, a=0..1),b=0..1), c = 0 .. 1),
   add(add(x[a,b],a=0..1),b=0..1)=1};
simplex[maximize](v, constraints, NONNEGATIVE); 

{v = 1/3, x[0,0] = 1/3, x[0,1] = 2/3, x[1,0] = 0, x[1,1] = 0}

Thus Player 1 should play your Action 2 with probability 1/3 and Action 4  with probability 2/3.  However, this does not maximize the expected payoff against any particular strategy of Player 2.  I don't understand what you're trying to do in your simulation.  What strategy is Player 2 supposed to be using?  If Player 2 is allowed to know Player 1's mixed strategy, you could simulate it as follows  This is using 1000 coin tosses for each data point.

X:= proc(pp)
      uses Statistics;
      local c,a,b, C, S, T, i;
      if pp <= 2/3 then c:= 0 else c:= 1 end if; # player 2's best strategy against pp
      S:= Sample(Bernoulli(pp),1000);  # player 1's moves on heads
      C:= Sample(Bernoulli(1/2),1000); # coin tosses
      T:= 0;
      for i from 1 to 1000 do
        if (C[i] = 0) and (S[i]=1) then T:= T-1 # player 1 passes on tails with prob pp
        elif c = 1 then T:= T+1  # player 2 passes
        else T:= T + 4*C[i]-2    # both play
        end if
     end do;
     T/1000
 end proc; 

randomize():
plot([seq([p,X(p)],p=0 .. 1, 0.01)],style=point); 

> plot3d(5,x=0..2,y=0..x^2,filled=true,axes=box,labels=[x,y,z],style=patchnogrid); 

Assuming your points represent a rectangular grid (m rows and n columns), you can use surfdata in the plots package.  There is no such thing as Matrixplot3d, as far as I know.
SurfacePlot produces three rows of points, one for each vector of data.  The x coordinates are 0 for the first data vector, 1/2 for the second and 1 for the third, the y coordinates go from 0 for the first data values to 1 for the last, and the z coordinates are the data values themselves.  Not being a statistician, I have no idea why anyone would want this sort of "surface plot".  But I find the help page description misleading - I too fell into the trap of thinking that this did what Christopher2222 was looking for.

That's not the error message I get.  I get 

Error, (in pdsolve/numeric/process_IBCs) initial/boundary conditions must depend upon exactly one of the independent variables: a(0,0) = 1

because your initial/boundary conditions are wrong: you could specify a(t,0) or a(0,y) but not just the single value a(0,0).

Isn't this the same system as in your previous question www.mapleprimes.com/forum/seperationvariables ?  See my reply there.

Are r, phi and theta all functions of t?  If you have three unknown functions, you need a system of three differential equations.

It is called a joke.

Your first problem is that H, Hy etc are expressions, not functions.  So it makes no sense to say e.g. H1 := diff(H(t,y), t).  It should be just H1 := diff(H, t).

It might help if you told us the actual code that produced this error. 

"Invalid terms in product" could mean that you are telling Maple to multiply something that isn't allowed to be multiplied.  For example, a string:

> "a" * b;
Error, invalid terms in product

Your PDE is second-order in x, so you need two boundary conditions.  Try it with, say, 
 

> BC:={T(0,t)=500*(1-t/1200)^10+1000,T(0.25,t) = 1500,T(x,0)=1500};

For stiff ODE systems, you might try dsolve(..., numeric, method=lsode).  See the help page

?dsolve,lsode

> X:= piecewise(t<2, 1.146967180*(t-1)^3-2.178399723*(2-t)^3+38.3612214*t+106.1461765, 
                     -.4672626810*(t-2)^3+1.146967180*(3-t)^3+45.2430246*t+92.3825701);
  Y:= piecewise(t<2, .3047483217*(t-1)^3-.5750322227*(2-t)^3+10.1163822*t+444.0230985,
                     -.1402738345*(t-2)^3+.3047483217*(3-t)^3+11.9448719*t+440.3661194);
  plot([X,Y, t = 1 .. 3]);