How to ad-hoc solve Einstein-equation for general metric numerically?

Question

Hey folks, feel free to move this post wherever you feel it fits (newbie)

I would like to solve some specific energy-momentum field configuration on a grid (say, 1+1 dim) and I wonder how this would be most efficiently done. Let us assume I have some initial data $T_{\mu \nu} (t,x)$ in some local coordinates $x$ at time $t = 0$, for instance some local energy density $\epsilon(x)$ and 4-speed $v_\mu(x)$, s.t.

$$T_{\mu\nu}(t=0, x) \equiv  \epsilon(x) v_\mu(x) v_\nu(x)$$

My first guess would be to start now by assuming that the system is slowly changing, i.e. $\partial_0 g_{\mu\nu} \approx 0$, assuming a minkowsky metric $\eta_{\mu\nu}$ with perturbations, calculating the Christoffel symbols, Curvature Tensor, etc. and finding the difference $\xi(x)$ at $t=0$:

$$\xi(x) = \left \|  8\pi T_{\mu\nu}(x) - R_{\mu\nu}(x) + g_{\mu\nu}(x)R(x)\right \|$$

Now I imagined using some crazy downhill algorithm to minimize this function, something like

$$g_{\mu\nu} \rightarrow g_{\mu\nu} - \frac{\partial \xi(x)}{\partial g_{\mu\nu}(x)}$$

Of course, this does not use the symmetries (such as Bianchi, etc.) and could fall into a local minimum (worse: Maybe the gauge condition $\nabla_\mu g_{\alpha \beta} = 0$ is not fullfilled?) However, let us imagine this worked, my next step would be to proceed to the next timestep $t \rightarrow t + \Delta t$ and use the geodetic equation $\nabla_\mu T^{\mu \nu} = 0$ to somehow calculate the new energy tensor $T^{\mu \nu}(t+\Delta t, x)$, i.e. assuming that locally we can do something badly forbidden like

$$\epsilon(x, t + \Delta t) = \frac{1}{V} \int d^d x' \cdot \epsilon(x',t) \delta(x,x' + \Delta t \cdot v(x')) \\ \tilde{v}(x, t + \Delta t) = \frac{1}{V} \int d^d x' \cdot v(x',t) \delta(x,x' + \Delta t \cdot v(x'))\\ v^\mu(x, t+ \Delta t) = \tilde{v}^\mu(x, t + \Delta t) - \Delta t \cdot \Gamma^\mu_{\alpha \beta} \tilde{v}^\alpha(x,t+\Delta t) \tilde{v}^\beta(x,t+\Delta t)$$

Where in the last equation I have used the geodetic equation for local speeds (?). After having calculated the new speeds and energy densities, the respective tensor follows as above and I can calculate the new metric. I think the main problem with the method here is the usage of the local frames $t$. Probably I have to change the coordinate system to the respective local frame to be able to make the last step and be able to use t?

Let us assume that this method has run through the whole world sheet, I would suggest just to start over and recalculate the metric, this time with changing derivates? Apparently the problem is to solve the partial differential equations with border condition $T_{\mu\nu}(t=0,x) = \text{const}$ and therefore the method described above is very naive. Maybe you have a better idea?