Solution (source code)

= Solution

For a <variation> $y+\varepsilon\eta$ with $\eta(0)=\eta(1)=0$, stationarity gives
$$
0=\int_0^1(F_y\eta+F_{y'}\eta')\,dx
=\int_0^1\left(F_y-\frac d{dx}F_{y'}\right)\eta\,dx.
$$
The boundary term vanishes. The <fundamental lemma of the calculus of variations> therefore gives the <Euler-Lagrange equation> $dF_{y'}/dx-F_y=0$.

Differentiating the proposed first integral gives
$$
\frac d{dx}(F-y'F_{y'})=F_x+y'\left(F_y-\frac d{dx}F_{y'}\right).
$$
Thus it is constant along an extremal when $F_x=0$, proving the <Beltrami identity>.

For the positive branch relevant to the endpoint data, $F=\sqrt{1+y'^2}/y$. The identity becomes
$$
\frac1{y\sqrt{1+y'^2}}=\frac1R,\qquad y^2(1+y'^2)=R^2.
$$
The <Euler-Lagrange equation> itself gives $yy''=-(1+y'^2)$, so $(y^2)''=-2$. Integrating, or solving the first integral, yields a circle $y^2+(x-c)^2=R^2$. The equal endpoint values force $c=1/2$ and $R^2=5/4$. Hence the unique positive smooth extremal is
$$
\boxed{y(x)=\sqrt{\frac54-\left(x-\frac12\right)^2}.}
$$
It stays positive on the entire interval and satisfies the second-order equation, including its turning point. Merely differentiating the first integral after dividing by $y'$ would need extra care there. Geometrically this is an arc of a <hyperbolic geodesic> in the <upper half-plane model>.