Boundary conditions at infinity for 1+2D wave equation in Mathematica 7

To solve a waves equation, I need to define some boundary conditions. The wave is propagating on an infinite plane, and it’s not a membrane fixed on some fixed support. I’m have difficulties in setting up the proper conditions with this working code:

Tmax = 20;

WaveEquation = 
  D[Phi[t, x, y], t, t] - D[Phi[t, x, y], x, x] - D[Phi[t, x, y], y, y] == 0;

PHI0[x_, y_] = Tanh[x] - 1.5Exp[-0.25(x^2 + y^2)]; (* Initial field *)

BoundaryConditions = { (* The issue is here ! *)
   Phi[t, -10, y] == PHI0[-10, y],
   Phi[t, 10, y] == PHI0[10, y],
   Phi[t, x, -10] == PHI0[x, -10],
   Phi[t, x, 10] == PHI0[x, 10]
   };

InitialConditions = {
   Phi[0, x, y] == PHI0[x, y],
   (D[Phi[t, x, y], t]/.t -> 0) == 0
   };

WaveSolution = NDSolve[
   Flatten@{WaveEquation, InitialConditions, BoundaryConditions},
    Phi,
   {t, 0, Tmax},
   {x, -10, 10},
   {y, -10, 10},
   PrecisionGoal -> 2
   ];

Manipulate[
 Plot3D[
  Evaluate[Phi[t, x, y]/.WaveSolution/.t -> time],
  {x, -10, 10},
  {y, -10, 10},
  PlotPoints -> {20, 20},
  PlotRange -> {{-10, 10}, {-10, 10}, {-3, 3}},
  Axes -> True,
  MeshFunctions -> (#3&),
  ColorFunction -> "Rainbow",
  ImageSize -> 700,
  Method -> {"RotationControl" -> "Globe"},
  SphericalRegion -> True
  ],
 {time, 0, Tmax, 0.01}
 ]

By running this code, we could clearly see the waves reflections on the bounding box. I want to get rid of that effect, since the wave is supposed to be propagating on an infinite plane. What can I do here? Take note that I’m working on Mathematica 7.0, so the option PeriodicBoundaryCondition isn’t recognized by my version, and it would be inapropriate anyway since the waves would get back to the origin after a while.

So how can I tell Mathematica that there isn’t any boundary to constraint the waves?

---

EDIT1: Is there a way to tell Mathematica that all waves should get truncated at the boundaries, without producing any reflections, creating something equivalent to a boundary at infinity? My problem mainly comes from the reflections at the borders {{-10, 10}, {-10, 10}} that I would like to prevent.

---

EDIT2: My question is similar to this one (but without an answer yet): Absorbing/Derivative Boundary Conditions

domain = {-10, 10}; tend = 20;

patt = x | y | t;
Table[{Numerator@
    Together@With[{theothervar = First@Complement[{x, y}, {var}]}, 
      var - PadeApproximant[#, {s, 0, {2, 0}}] & /@ 
        Simplify[var /. Solve[t^2 == x^2 + y^2, var] /. theothervar -> s t, t > 0] /. 
       s -> theothervar/t], var == Reverse@domain // Thread}[Transpose], {var, {x, y}}]

bc = % /. {(var1 : patt) (var2 : patt) :> D[u[t, x, y], var1, var2], 
           (var : patt)^n_ :> D[u[t, x, y], {var, n}]}

With[{u = u[t, x, y]}, 
 ic = {D[u, t] == 0, u == Tanh[x] - 1.5 Exp[-0.25 (x^2 + y^2)]} /. t -> 0;
 eq = Piecewise[Flatten[bc, 1], D[u, t, t] - (D[u, x, x] + D[u, y, y])] == 0]

points = 50; difforder = 2;
grid = Rescale[Range@points, {1, points}, domain];

(*Definition of pdetoode isn't included in this post,
  please find it in the link above.*)
ptoofunc = pdetoode[u[t, x, y], t, {grid, grid}, difforder];
ode = ptoofunc@Simplify`PWToUnitStep@eq;
odeic = ptoofunc@ic;

sollst = Partition[
   NDSolve[{ode, odeic}, Flatten@Outer[u, grid, grid], {t, 0, tend}][[1, All, -1]], 
   points];

solfunc = rebuild[sollst, {grid, grid}]

Table[Plot3D[solfunc[t, x, y], {x, ##}, {y, ##}, PlotRange -> {-1.5, 1.5}] & @@ 
   domain, {t, 0, 20, 1/2}] // ListAnimate