Jacobi equation for the geodesic deviation in the weak field limit

The geodesic deviation equation can be written in the following form
$$
nabla_U^2 xi = R (U, xi) U
tag{1}
$$
where $R$ is the Ricci tensor. It can also be written component-wise using the Riemann tensor
$$
(nabla_U^2 xi)^alpha = R^alpha_{;, beta mu nu} : U^beta U^mu xi^nu
tag{2}
$$

I’m interested in this equation on a Riemannian manifold equipped with the following metric tensor
$$
g = -(1+2phi(x))mathrm{d} t otimes mathrm{d} t + (1-2phi (x))left(mathrm{d} x otimes mathrm{d} x + mathrm{d} y otimes mathrm{d} y + mathrm{d} z otimes mathrm{d} z right)
tag{3}
$$
with the standard torsion-free and metric-preserving connection $nabla$. I’m only interested in the spatial part of the Jacobi equation, that starts from $U = e_0$ i.e.
$$
(nabla_t^2 xi)^i = R^i_{;, 0 0 nu} ,xi^nu
tag{4}
$$

In textbooks it can be found that in the weak-field limit (keeping only the linear terms involving $phi (x)$) this reduces to
$$
frac{mathrm{d}^2 xi^i}{mathrm{d} t^2} = – phi_{,ij} , xi^j
tag{5}
$$
where comma indicates partial derivatives.

However, when I try this for $g$ above, I get a different result. First, the left-hand side of (4) is
$$
nabla_t xi = left( xi^mu_{,t} + Gamma^mu_{;; nu t} , xi^nu right) e_mu
$$
$$
nabla^2_t xi = left( xi^mu_{,t} + Gamma^mu_{;; nu t} , xi^nu right)_{,t} e_mu + left( xi^mu_{,t} + Gamma^mu_{;; nu t} , xi^nu right) Gamma^lambda_{;; mu t} e_lambda
$$

The Christoffel symbols are in general
$$
Gamma^alpha_{;; mu nu} = frac{1}{2} g^{alpha lambda} left( g_{lambda mu, nu} + g_{lambda nu, mu} – g_{mu nu, lambda} right)
$$
so in our case
$$
Gamma^0_{;; mu nu} = begin{pmatrix}
phi_{,0} & phi_{,1} & phi_{,2} & phi_{,3}
phi_{,1} & – phi_{,0} & 0 & 0
phi_{,2} & 0 & – phi_{,0} & 0
phi_{,3} & 0 & 0 & – phi_{,0}
end{pmatrix}
$$
$$
Gamma^1_{;; mu nu} = begin{pmatrix}
phi_{,1} & – phi_{,0} & 0 & 0
– phi_{,0} & – phi_{,1} & – phi_{,2} & – phi_{,3}
0 & – phi_{,2} & phi_{,1} & 0
0 & – phi_{,3} & 0 & phi_{,1}
end{pmatrix}
$$
$$
Gamma^2_{;; mu nu} = begin{pmatrix}
phi_{,2} & 0 & – phi_{,0} & 0
0 & phi_{,2} & – phi_{,1} & 0
– phi_{,0} & – phi_{,1} & – phi_{,2} & – phi_{,3}
0 & 0 & – phi_{,3} & phi_{,2}
end{pmatrix}
$$
$$
Gamma^3_{;; mu nu} = begin{pmatrix}
phi_{,3} & 0 & 0 & – phi_{,0}
0 & phi_{,3} & 0 & – phi_{,1}
0 & 0 & phi_{,3} & – phi_{,2}
– phi_{,0} & – phi_{,1} & – phi_{,2} & – phi_{,3}
end{pmatrix}
$$
which, for $(nabla^2_t xi)^i$ gives (to the linear order in $phi$)
$$
xi^i_{;;, 00} + phi_{, i0} xi^0 – phi_{,00} xi^i + 2 phi_{,i} xi^0_{;;, 0} – 2 phi_{,0} xi^i_{;;,0}
$$

The right-hand side needs $R^i_{;; 00 mu}$. First, due to the symmetries, $R^i_{;;000} = 0$, so we only need
$$
R^i_{;; 00j} = – phi_{, ij} – phi_{,00} , delta_{ij}
$$

Now put it all together
$$
frac{partial^2 xi^i}{partial t^2} + phi_{, i0} xi^0 – phi_{,00} xi^i + 2 phi_{,i} xi^0_{;;, 0} – 2 phi_{,0} xi^i_{;;,0} = – phi_{,ij} xi^j – phi_{,00} xi^i
$$

Cancel out common terms
$$
frac{partial^2 xi^i}{partial t^2} + phi_{, i0} xi^0 + 2 phi_{,i} xi^0_{;;, 0} – 2 phi_{,0} xi^i_{;;,0} = – phi_{,ij} xi^j
$$

There are problems with this expression. It does not match what the textbooks say.

Even if it did, how can I get "$frac{mathrm{d}^2 xi^i}{mathrm{d} t^2}$" instead of the partial derivative? Is $frac{mathrm{d}}{mathrm{d} t}$ to be interpreted as $U^mu partial_mu$? In that case if $U = e_0$, we would have $frac{mathrm{d}}{mathrm{d} t} = frac{partial}{partial t}$, but probably only at the initial point of our geodesic…or do we reparametrize the geodesic in terms of the time $t$ and then derivative w.r.t $t$ is actually derivative w.r.t. the curve parameter, so $nabla_U^2 = frac{mathrm{d}^2}{mathrm{d} lambda^2} = frac{mathrm{d}^2}{mathrm{d} t^2}$? That would certainly interpret the left-hand side of (4), but we would still have that weird term $-phi_{,00} xi^i$ on the right-hand side (unless we assume that $phi$ does not depend on time, on the top of all that). But even if we do it, then $U$ is no longer in just zeroth direction so the right-hand side will feature two $U$ terms (even if we begin with $U = e_0$, can it stay that way for the entire geodesic? or do we assume small $v$, therefore the zeroth component dominates and since $U cdot U = -1$, then $U^0 = 1$?)

There’s a bunch of other terms that don’t seem to arise in textbooks, but I never saw the expression (5) derived, only stated.

Is there a conceptual mistake in my thinking, or I miscalculated something?