Notation:
A^a_,b = ∂A^a/∂x^b and A_i,j = ∂A_i/∂x^j
A^a_;b = ∇_bA^a

The idea here is to write the simplest program to carry out the computation explicit in these equations:
∂_tg_ab = −2αK_ab + ∇_aβ_b + ∇_bβ_a
∂_tK_ab = −∇_a∇_bα + α(R_ab + K_abK − 2K_acK^c_b) + β^c∇_cK_ab + K_ca∇_bβ^c + K_cb∇_aβ^c
By simplest I presume to minimize the effort to understand the program and see that it is indeed doing the intended calculation.

We must compute the Christoffel symbols:
Γⁱ_kl = (1/2)g^im(g_mk,l + g_ml,k − g_kl,m)
and for that we need g^im — the inverse of the metric tensor.

With these values the covariant derivative is:
∇_aβ_b = β_b;a = ∂β_b/∂x^a − Γ^c_abβ_c
along with the other familiar equations.

Taking simple derivatives of field values with respect to spacial coordinates is fraught with centering issues. I plan to use the simplest centering, j, k, l < n} schemes and document how each value is centered within the computational zone.

A significant check of code such as this is to use Schwarzschild or even Kerr metrics as initial conditions.

For familiarity I assume calculations on a mesh [0:n]⁴ = {(i, j, k, l) | 0 ≤ i, j, k, l < n}. Posit: g_ab values kept for lattice points. Tentatively we keep K_ab there too.