1+1 Numerical Relativity Simulation

I am trying to teach myself numerical relativity by starting with the simplest non-trivial scenario: a 1+1 vacuum spacetime with a non-trivial initial slice. Essentially, I am following this paper and trying to reproduce the results for the flat 2d (1+1) spacetime. The problem is also discussed on page 364 of "Introduction to 3+1 relativity".

Since there is only one spatial dimension, let $$g\equiv g_{xx}$$ and $$K\equiv K_{xx}$$. Also take the shift $$\beta^{i}$$ to be zero. The system of PDEs to be evolved is made first order by defining $$D_{\alpha}\equiv\partial_{x}\ln \alpha$$, $$D_{g}\equiv\partial_{x}\ln g$$, and $$\tilde{K}\equiv\sqrt{g}K$$. Thus the system consists of the five evolving fields:

$$\partial_{t}\alpha = -\alpha^{2}f\frac{\tilde{K}}{\sqrt{g}} \\ \partial_{t}g = -2\alpha\sqrt{g}\tilde{K} \\ \partial_{t}D_{\alpha} = -\partial_{x}\left(\alpha f\frac{\tilde{K}}{\sqrt{g}}\right) \\ \partial_{t}D_{g} = -\partial_{x}\left(2\alpha\frac{\tilde{K}}{\sqrt{g}}\right) \\ \partial_{t}\tilde{K} = -\partial_{x}\left(\frac{\alpha D_{\alpha}}{\sqrt{g}}\right)$$ Here, I am using $$f=1$$ for the harmonic slicing condition.

The spacetime is a vacuum, but the problem is to study the gauge dynamics of using a non-trivial initial slice. Such a slice may be defined in the Minkowski coordinates:

$$t_{M}=h(x_{M})$$

Here, $$h$$ is chosen to be a Gaussian. The lapse $$\alpha$$ is taken to be initially 1 everywhere.

Therefore, the initial value problem is $$\alpha(0,x) = 1 \\ g(0,x) = 1 - h'^{2}\\ D_{\alpha}(0,x) = 0 \\ D_{g}(0,x) = \frac{2h'h''}{g} \\ K(0,x) = -\frac{h''}{g} \\$$

From this I discretize the system and advance all the fields simultaneously in an FTCS scheme. (I know it's unstable but I want to get it working before moving on to a more advanced scheme.)

The results are shown in the paper. Basically what is supposed to happen is that in every field, there develop two wave pulses propagating in either direction. They should travel at speed $$\sqrt{f}=1$$.

However, in my case, I get a sort of waveform that appears immediately but does not propagate, instead it just increases in amplitude. I am fairly sure there are no errors in my code, so I believe I am missing something conceptually. I have the right initial conditions which are also shown in the paper.

What puzzles me is that the author remarks that there should be "3 fields that propagate along the time lines (speed zero)". Two of these mentioned are $$\alpha$$ and $$g$$. Does this mean there is some sort of coordinate transformation that I need to make before trying to visualize the data?

Does anyone know of any explicit 1+1 numerical relativity routines that I could consult? I would like to see the actual code.

Well this is embarrassing. It turns out there was a very subtle typo in my code. So there is no real physical problem at all, and the approach described is accurate and there is no need for any coordinate transformations. However, I still found that I could not reproduce the same results in the paper. My propagating solutions were traveling much slower than what is expected, even when I used similar initial parameters for $$h(x)$$. I also discovered that the propagation speed in the lapse changed depending on my choice $$\Delta t$$ and $$\Delta x$$. Therefore, the problem stems from the first-order FTCS scheme chosen; it is not accurate enough.
For $$h(x) = \mathrm{e}^{-\frac{x^{2}}{\sigma^{2}}}\\ \sigma = 10.0$$ And using the same discretization used in the paper, $$\Delta t=0.125$$ and $$\Delta x=.25$$, I find:
From this it is clear that the pulses do not travel at $$\sqrt{f}=1$$, instead it would appear to be closer to $$4$$. Note that since the lapse is simply a gauge function there would be no problem with the pulses traveling at light speed.