Where does Jacobi's accessory equation come from? I'm reading Charles Fox's An Introduction to the Calculus of Variations and in section 2.4 he just suddenly introduces Jacobi's accessory equation and I don't understand where it's coming from.
Jacobi's accessory equation (which oddly doesn't have a Wikipedia page) is $$\left[\frac{\partial^2 F(x,s,s')}{\partial s^2}-\frac d{dx}\frac{\partial^2 F(x,s,s')}{\partial s\partial s'}\right]u-\frac d{dx}\left(\frac{\partial^2 F(x,s,s')}{\partial s'^2}\frac {du}{dx}\right)=0$$ and has something to do with the second variation of $\int_a^b F(x,y,y')dx$ evaluated near a stationary path $y=s(x)$.
Could someone either derive (or give the main idea of the derivation) the equation or suggest a good source that does?
P.S. I checked in Gelfand and Fomin's book which I'm not reading but have laying around. The problem with their derivation is that it's in the middle of their book (as opposed to near the beginning of Fox's) and thus seems to make use of stuff I haven't gotten to yet in Fox's book.

Here's the context in Fox's book:
Note that for brevity the notation $F_{00} := \frac{\partial^2 F}{\partial s^2}$, $F_{01} := \frac{\partial^2 F}{\partial s\partial s'}$, $F_{11} := \frac{\partial^2 F}{\partial s'^2}$ is used in the following.
We just started looking at the second variation and trying to derive the conditions for extremizing the functional subject to weak variation.  In the previous section we derived that if $t(a)=t(b)=0$, then $$\int_a^b \left(t^2F_{00} +2tt'F_{01}+t'^2F_{11}\right)dx = \int_a^b\left[t^2F_{00}-t^2\frac d{dx}(F_{01})-t\frac d{dx}(t'F_{11})\right]dx$$
Then this section starts off with:

On solving the [Euler-Lagrange equation], the equation of the extremal $y=s(x)$ which passes through the given points $A$ and $B$ can be determined.
Thus the quantities $F_{00}$, $F_{01}$, $F_{11}$, and $\frac d{dx} F_{01}$ can all be expressed in terms of $x$ and the differential equation $$\left[F_{00}-\frac d{dx}F_{01}\right]u-\frac d{dx}\left(F_{11}\frac {du}{dx}\right)=0 \tag{1}$$
can then be solved for $u$ as a function of $x$.  This is an ordinary linear differential equation of the second order and is known as the subsidiary or Jacobi's equation or, more frequently, as the accessory equation.
On taking $x$ to be independent and $t(=t(x))$ the dependent variable in the integral $I_2 [= \int_a^b \left(t^2F_{00} +2tt'F_{01}+t'^2F_{11}\right)dx]$, it is easily seen that $(1)$ is the [Euler-Lagrange equation] for minimizing $I_2$ with $t$ replaced by $u$.

Note that other than that last line he doesn't explain what the function $u$ is.
 A: There are at least three ways of getting to this equation: the first is to add zero to the second variation of the functional, in the form
$$ 0 = \int_a^b (wt^2)' \, dx = \int_a^b (2wtt'+ w't^2) \, dx, $$
since $t$ satisfies the boundary conditions $t(a)=t(b)=0$. Then the second variation becomes
$$ \int_a^b \left( \left( F_{00}-\frac{dF_{01}}{dx} \right)t^2 + 2wtt'+ (w'+F_{11})t'^2 \right) dx. $$
Write $p=F_{00}-\frac{dF_{01}}{dx}$, $q=F_{11}$ for brevity. Completing the square on the $t^2$ term gives
$$ \int_a^b \left( p\left( t+\frac{wt'}{p} \right)^2 + \left( w'-\frac{w^2}{p}+q \right)t'^2 \right) dx $$
It is easy to see that the first term is positive for any $w \neq 0$ (just solve the differential equation $t+\frac{wt'}{p}=0$) Therefore, the second variation will be positive provided that the second term is zero, which amounts to solving the equation
$$ w'-\frac{w^2}{p}+q = 0 $$
This is a Riccati equation, and we can make it into a linear second-order equation by the substitution $w=-pu'/u$, and then we find
$$ -(pu')' + qu = 0, $$
which you'll find is the Jacobi accessory equation; there are no boundary conditions on $u$, but $w$ only makes sense if $u$ has no zero in $[a,b]$.

(Alternatively, the Euler–Lagrange equation for the second variation being zero is the Jacobi accessory equation, but this approach doesn't explain the roots or the boundary conditions.)

The more sophisticated version derives the Jacobi accessory equation from the difference between two adjacent solutions to the Euler–Lagrange equations: suppose that $y$ and $y+hv$ solve the Euler–Lagrange equations, and satisfy the boundary conditions up to order $h$. Then we look for the equation satisfied by $v$. Since they both satisfy the equations, we have
$$ \begin{align}
0 &= \left( \frac{\partial L[x,y+hv,y'+hv']}{\partial y}-\frac{d}{dx}\frac{\partial L[x,y+hv,y'+hv']}{\partial y'} \right) - \left( \frac{\partial L[x,y,y']}{\partial y}-\frac{d}{dx}\frac{\partial L[x,y,y']}{\partial y'} \right) \\
&= h\left( v\frac{\partial^2 L}{\partial y^2} + v'\frac{\partial^2 L}{\partial y\partial y'}- \frac{d}{dx} \left( v \frac{\partial^2 L}{\partial y \partial y'} + v' \frac{\partial^2 L}{\partial y'^2} \right) \right) +o(h) \\
&= h\left( v\left(\frac{\partial^2 L}{\partial y^2}- \frac{d}{dx} \frac{\partial^2 L}{\partial y \partial y'} \right) - \frac{d}{dx} \left( v' \frac{\partial^2 L}{\partial y'^2} \right) \right) + o(h)
\end{align} $$
by cancelling one of the terms in the differentiation of the product. Dividing by $h$ and taking the limit as $h \to 0$ gives the Jacobi equation.
Here, $v$ has no boundary conditions specified, and $v$ having a zero $a$ is interpreted as the solutions all passing through the same point at $a$, which is called a conjugate point.
