Assume the initial condition is that a point mass starts at height $y_0$. After descending to height $y < y_0$, we know that its speed will be $v = \sqrt{2mg(y_0 - y)}$. Thus, the displacement element can be written as $$ds^2 = dx^2 + dy^2 = v^2 dt^2,$$ so that we have $$dt^2 = \frac{1}{2mg(y_0 - y)} (dx^2 + dy^2).$$ We are trying to minimize the functional $\int \sqrt{dt^2}$. This looks just like the geodesic problem (from general relativity) where the metric is $$g_{\mu \nu} = \frac{1}{2mg}\begin{pmatrix} \frac{1}{y_0 - y} & 0 \\ 0 & \frac{1}{y_0 - y}\end{pmatrix}.$$ The geodesic equation reads, $$\frac{d^2 x^\mu}{dt^2} + \frac{1}{2}\Gamma^{\mu}_{\alpha\beta} \frac{dx^\alpha}{dt} \frac{dx^\beta}{dt} = 0.$$
I believe that with this metric, we have the Christoffel symbols $$\Gamma^x_{xy} = \Gamma^x_{yx} = \frac{1}{2(y_0-y)}$$ $$\Gamma^y_{yy} = \frac{1}{2(y_0-y)}$$ $$\Gamma^y_{xx} = \frac{-1}{2(y_0-y)}$$ (note the negative sign)
and all else are zero.
This would give the equations $$\frac{d^2 x}{dt^2} + \frac{1}{2(y_0-y)}\frac{dx}{dt}\frac{dy}{dt} = 0$$
$$\frac{d^2 y}{dt^2} + \frac{1}{2(y_0-y)}\left(\frac{dy}{dt}\right)^2 - \frac{1}{2(y_0-y)}\left(\frac{dx}{dt}\right)^2= 0.$$
Is this all correct so far, and can this be continued into a useful solution to the problem? I am also concerned that I might not have the right derivative in the geodesic equation - I think it should be derivatives with respect to $t$ but I could imagine being wrong and it should be with $s$, related by $$ds = vdt \implies \frac{d}{dt} = v \frac{d}{ds}.$$
We know the solution for the brachistochrone problem is a cycloid and we have a simple parametric form for those; we can subtitute into the geodesic equation you wrote down to see if they are satisfied by the cycloid. If so, that establishes that the approach is valid.
To simplify the calculation, let's take $y_0 =0$, reverse gravity so that it is uniform but pointing upward, set $2g=1$, so the motion occur in the upper half plane (moving upward from the x-axis). One other minor detail: in deriving the equation, equating the kinetic and potential energy, the mass terms cancel and so with these conventions $v = \sqrt{y}$. Then the metric you are considering then takes the form $\frac{dx^2+dy^2}{y}$. To compute the equations of the geodesics we use the following formulas: \begin{align} [ij,k] &= \frac{1}{2}(\partial_i g_{jk} +\partial_j g_{ik} - \partial_kg_{ij}) \\ \Gamma^i_{jk} &= g^{il}[jk,l] \end{align} The equation for the geodesics are then: \begin{equation} \frac{d^2x^i}{dt^2} + \Gamma^i_{jk} \frac{dx^j}{dt}\frac{dx^k}{dt}= 0 \end{equation}
The parametric equations for the cycloid that solves the brachisochrone problem (with $r>0$ constant) are: \begin{align} x &=r(t-\sin(t)) & \dot{x} &= r(1-\cos(t)) & \ddot{x} &= r\sin(t)\\ y &=r(1-\cos(t)) & \dot{y} &= r\sin(t) & \ddot{y} &=r\cos(t) \end{align} Note that $\dot{x} = y$.
Specializing with coordinates $\{x,y\}$: \begin{align} \ddot{x} + \Gamma^x_{xx}\dot{x}^2 +2\Gamma^x_{xy}\dot{x}\dot{y} +\Gamma^x_{yy}\dot{y}^2 &= 0 \\ \ddot{y} + \Gamma^y_{xx}\dot{x}^2 +2\Gamma^y_{xy}\dot{x}\dot{y} +\Gamma^y_{yy}\dot{y}^2 &= 0 \end{align}
and for $g^{-1}=\begin{pmatrix} g^{xx} & g^{xy} \\ g^{yx} & g^{yy}\end{pmatrix} = \begin{pmatrix} y & 0 \\ 0 & y \end{pmatrix}$
We get : \begin{align*} \Gamma^x_{xx} &= g^{xx}[xx,x] = 0 \\ \Gamma^x_{xy} &= g^{xx}[xy,x] = -1/2y \\ \Gamma^x_{yy} &= g^{xx}[yy,x] = 0 \end{align*}
Similarly we compute: \begin{align*} \Gamma^y_{xx} &= g^{yy}[xx,y] = 1/2y \\ \Gamma^y_{xy} &= g^{yy}[xy,y] = 0 \\ \Gamma^y_{yy} &= g^{yy}[yy,y] = -1/2y \end{align*}
(note the signs on the $\Gamma^y_{xx} \text{ and } \Gamma^y_{yy}$ are opposite yours)
This yields the two equations: \begin{align} \ddot{x} - \frac{1}{y}\dot{x}\dot{y} &= 0 & (1)\\ \ddot{y} + \frac{1}{2y} \dot{x}^2 - \frac{1}{2y}\dot{y}^2 &= 0 & (2) \end{align}
Eq (1) is easily seen to be true using the fact that $\dot{x} = y$ and (2) is also verified with some algebra. So I conclude that the approach works!