The numerical solution of the initial value problems (2.1) - (2.2) and (2.5) can efficiently and conveniently be done by a finite difference method.
In view of our stability result (see previous section) this suggestions itsself, but it has been used already in [9].
We simply use the usual five point discretization on the grid
where 
and 
and 
. Then, the finite difference approximation to (2.1) - (2.2) reads (the superscript j is omitted)
The boundary conditions in (3.1) assume that the field can be measured everywhere on each of the faces of 
.
The values of 
for 
have to be computed from the field in 
by numerical differentiation. (3.1) is solved in a recursive way, i.e. if w in known on the levels 
, 
we compute it for 
by (3.1). In order to preserve stability we have to filter out the frequencies greater than 
at each stage of this process. This can be done by doing a 2D FFT on each matrix 
as soon as it is computed, zeroing all the entries with 
, followed by a 2D inverse FFT. The total operator count for solving (2.1) - (2.2) once is O 
. For details see [16]. (2.5) is solved exactly in the same way.
