So using Fermat's principle, the light path should be the path which takes the least amount of time. We need to have a picture, so that we know exactly what it is we are doing (or at least so you know what I'm doing). Now, I'm going to write out the equations and solve them (because I don't think your's is correct). What I'm going to find is that it's difficult to solve, but changing the function I'm solving for will make it easier (but I'm still going to write out the "hard" one for those who don't see why you might want to try this).
$$
|d\vec{s}| = |\vec{v}|dt \rightarrow \sqrt{dx^2 + dy^2} = |\vec{v}|dt \\
dy = f'dx \rightarrow |d\vec{s}| = \sqrt{1 + f'^2}dx\\
t_\text{total} = \int dt = \int \frac{\sqrt{1 + f'^2}}{|\vec{v}|}dx \\
t_\text{total} = \frac{1}{c_0}\int \sqrt{1 + f'^2}\sqrt{1 - kf}dx
$$
Now find the Euler-Lagrange equation:
\begin{align}
\mathcal{L} =& \sqrt{\left(1 + f'^2\right)\left(1 - kf\right)} \\
\frac{\partial \mathcal{L}}{\partial f} = &-\frac{k}{2}\sqrt{\frac{1 + f'^2}{1 - kf}} \\
\frac{\partial \mathcal{L}}{\partial f'} = &f'\sqrt{\frac{1 - kf}{1 + f'^2}} \\
\frac{d}{dx}\frac{\partial \mathcal{L}}{\partial f'}=& f''\sqrt{\frac{1 - kf}{1 + f'^2}} - \frac{kf'^2}{2\sqrt{\left(1 - kf\right)\left(1 + f'^2\right)}} - \frac{f'^2f''}{\sqrt{\left(1 + f'^2\right)^3\left(1 - kf\right)}}
\end{align}
So this is huge mess. Instead of trying to find $y = f(x)$, let's try to find $x = f(y)$. This doesn't change the $|d\vec{s}|$ except that now it's $\sqrt{1 + f'^2}dy$. Now the new Lagrangian is:
\begin{align}
\mathcal{L} =& \sqrt{1 + f'^2}\sqrt{1 - ky} \\
\frac{\partial \mathcal{L}}{\partial f} =& 0\\
\frac{\partial \mathcal{L}}{\partial f'} =& f'\sqrt{\frac{1 - ky}{1 + f'^2}}
\end{align}
Now, since this gives the following equation:
$$
\frac{d}{dy}\left(f'\sqrt{\frac{1 - ky}{1 + f'^2}}\right) = 0
$$
We can just set that partial to some constant:
$$
f'\sqrt{\frac{1 - ky}{1 + f'^2}} = \alpha \\
f'^2\left(1 - ky\right) = \alpha^2\left(1 + f'^2\right) \\
f'^2\left(1 - ky - \alpha^2\right) = \alpha^2 \\
f' = \frac{\alpha}{\sqrt{1 - ky - \alpha^2 }}
$$
Now integrate both sides:
$$
f = -\frac{2}{k}\alpha\sqrt{1 - ky - \alpha^2} + \beta
$$
Go ahead and sub back in $f = x$ and solve for $y$:
$$
x - \beta = -\frac{2}{k}\alpha\sqrt{1 - ky - \alpha^2} \\
\frac{k^2(x - \beta)^2}{4\alpha^2} = 1 - ky - \alpha^2 \\
y = \frac{1 - \alpha^2 - \frac{k^2(x - \beta)^2}{4\alpha^2}}{k} \\
y = \frac{4\alpha^2(1 - \alpha^2) - k^2(x - \beta)^2}{4k\alpha^2} \\
y = \frac{4\alpha(1 - \alpha) - k^2(x - \beta)^2}{4k\alpha}
$$
In the last step, I changed from $\alpha^2$ to just $\alpha$--so those last two equations aren't strictly equal, but since $\alpha$ is a constant, I can make it whatever I want (although this should mean $\alpha \geq 0$). Now you need to find $\alpha$ and $\beta$ for your boundary conditions:
$$
f(-x_0) = 0 = \frac{4\alpha(1 - \alpha) - k^2(-x_0 - \beta)^2}{4k\alpha} = \frac{4\alpha(1 - \alpha) - k^2(x_0 + \beta)^2}{4k\alpha} \\
f(x_0) = 0 = \frac{4\alpha(1 - \alpha) - k^2(x_0 - \beta)^2}{4k\alpha} \\
4\alpha(1 - \alpha) - k^2(x_0 + \beta)^2 = 0 \\
4\alpha(1 - \alpha) - k^2(x_0 - \beta)^2 = 0 \\
(x_0 + \beta)^2 = (x_0 - \beta)^2 \rightarrow \beta = 0 \\
4\alpha(1 - \alpha) = k^2x_0^2 \\
f(x) = \frac{k^2}{\zeta}\left(x_0^2 - x^2\right)
$$
Here $\zeta = 4k\alpha$ for whatever value $\alpha$ attains by solving the equation above it. Hopefully it's obvious that f(x) attains a maximum when $x = 0$ (the vertex of the parabola centered around $x = 0$) and that this maximum is $y_0 = \frac{k^2x_0^2}{\zeta} = \frac{k^2x_0^2}{4k\alpha}$. Now substitute $\alpha = \frac{kx_0^2}{4y_0}$ into the equation we never solved:
$$
4\alpha(1 - \alpha) = k^2x_0^2 \rightarrow 4\frac{kx_0^2}{4y_0}\left(1 - \frac{kx_0^2}{4y_0}\right) = k^2x_0^2 \\
\frac{k}{y_0}\cdot\frac{4y_0 - kx_0^2}{4y_0} = k^2 \\
4y_0 -kx_0^2 = 4ky_0^2 \\
kx_0^2 = 4y_0 - 4ky_0^2 = 4y_0(1 - 4ky_0)
$$
Now just multiply both sides by $k$:
$$
k^2x_0^2 = 4ky_0\left(1 - 4ky_0\right) \rightarrow (kx_0)^2 = 4ky_0\left(1 - ky_0\right) \text{, q.e.d}
$$
