I'm trying to solve a very typical wave equation as a test problem. The codes are given below,

L = 10;
\[Lambda]  = 2^0 L;
c = 1;
\[Omega] = 2 \[Pi] (c/L) (1/100);
tMax = 100 L/c;

Equation =  D[\[Phi][x, t], t, t] == c^2 D[\[Phi][x, t], x, x];
IConds = {\[Phi][x, 0] == 0, Derivative[0, 1][\[Phi]][x, 0] == 0};
BConds = {Derivative[1, 0][\[Phi]][0, t] == Tanh[t] Sin[\[Omega] t], 
   Derivative[1, 0][\[Phi]][L, t] == 0};

Sol = NDSolve[
  Join[{Equation}, BConds, IConds], \[Phi], {x, 0, L}, {t, 0, 
   tMax}, Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid",  
      "MinPoints" -> 100}}]

The problem can be solved for up to t ~ 10 L/c; for t > 20 L/c, the solution grows due to instabilities. I'm not sure about the discretization procedures that Mathematica employed in solving hyperbolic PDEs; it seems that if discretization had been done correctly the solution would not have blown up. For example, in one type of discretization scheme the Courant condition c/(Delta x/Delta t) << 1 has to be satisfied, where Delta x is the gridsize in the spatial direction, Delta t the gridsize in the temporal direction, and c the speed of wave.

My question is just simply, is there a way for Mathematica to integrate a wave-like equation in the temporal direction stably, over a long time t, that is, t = n T, where n is an integer and T = L/c?