[proofplan]
We compute the first and second derivatives of the WKB ansatz $\psi_\hbar = ae^{iS/\hbar}$ exactly. Substituting the second derivative into the stationary Schrödinger equation factors out the nonzero exponential $e^{iS/\hbar}$ and leaves an expansion in powers of $\hbar$. The coefficients of $\hbar^0$ and $\hbar^1$ are then identified directly, giving the Hamilton-Jacobi and transport equations. Finally, the transport equation is rewritten as a conservation law for $a^2S'$.
[/proofplan]
[step:Differentiate the WKB ansatz twice]
Fix $\hbar > 0$. Define the phase factor $q_\hbar: I \to \mathbb{C}$ by
\begin{align*}
q_\hbar(x) = e^{iS(x)/\hbar}.
\end{align*}
Since $S \in C^2(I;\mathbb{R})$, the chain rule gives
\begin{align*}
q_\hbar'(x) = \frac{i}{\hbar}S'(x)q_\hbar(x).
\end{align*}
Since $\psi_\hbar(x)=a(x)q_\hbar(x)$ and $a \in C^2(I;\mathbb{C})$, the product rule gives
\begin{align*}
\psi_\hbar'(x) = \left(a'(x) + \frac{i}{\hbar}a(x)S'(x)\right)q_\hbar(x).
\end{align*}
Differentiating once more and again using the product and chain rules,
\begin{align*}
\psi_\hbar''(x) = \left(a''(x) + \frac{2i}{\hbar}a'(x)S'(x) + \frac{i}{\hbar}a(x)S''(x) - \frac{1}{\hbar^2}a(x)(S'(x))^2\right)q_\hbar(x).
\end{align*}
[guided]
Fix $\hbar > 0$. The only $\hbar$-dependence in the ansatz comes from the oscillatory exponential, so we isolate it by defining $q_\hbar: I \to \mathbb{C}$ by
\begin{align*}
q_\hbar(x) = e^{iS(x)/\hbar}.
\end{align*}
Because $S \in C^2(I;\mathbb{R})$, the chain rule applies and gives
\begin{align*}
q_\hbar'(x) = \frac{i}{\hbar}S'(x)q_\hbar(x).
\end{align*}
Now use $\psi_\hbar(x)=a(x)q_\hbar(x)$. Since $a \in C^2(I;\mathbb{C})$, the product rule gives
\begin{align*}
\psi_\hbar'(x) = a'(x)q_\hbar(x) + a(x)q_\hbar'(x).
\end{align*}
Substituting the formula for $q_\hbar'$ yields
\begin{align*}
\psi_\hbar'(x) = \left(a'(x) + \frac{i}{\hbar}a(x)S'(x)\right)q_\hbar(x).
\end{align*}
We differentiate this expression once more. The derivative of the bracketed factor is
\begin{align*}
a''(x) + \frac{i}{\hbar}a'(x)S'(x) + \frac{i}{\hbar}a(x)S''(x).
\end{align*}
The derivative of $q_\hbar$ contributes an additional factor $\frac{i}{\hbar}S'(x)q_\hbar(x)$. Therefore
\begin{align*}
\psi_\hbar''(x) = \left(a''(x) + \frac{i}{\hbar}a'(x)S'(x) + \frac{i}{\hbar}a(x)S''(x)\right)q_\hbar(x) + \left(a'(x) + \frac{i}{\hbar}a(x)S'(x)\right)\frac{i}{\hbar}S'(x)q_\hbar(x).
\end{align*}
Collecting the two $a'S'$ terms and using $i^2=-1$ gives
\begin{align*}
\psi_\hbar''(x) = \left(a''(x) + \frac{2i}{\hbar}a'(x)S'(x) + \frac{i}{\hbar}a(x)S''(x) - \frac{1}{\hbar^2}a(x)(S'(x))^2\right)q_\hbar(x).
\end{align*}
This exact second-derivative formula is the whole source of the WKB hierarchy.
[/guided]
[/step]
[step:Substitute into the Schrödinger equation and collect formal powers of $\hbar$]
For each $x \in I$, substitute the expression for $\psi_\hbar''(x)$ into the residual map $R_\hbar:I\to\mathbb{C}$ defined by
\begin{align*}
R_\hbar(x) = -\frac{\hbar^2}{2m}\psi_\hbar''(x) + V(x)\psi_\hbar(x) - E\psi_\hbar(x).
\end{align*}
Since $q_\hbar(x)=e^{iS(x)/\hbar}$ is nonzero, the residual equals $q_\hbar(x)$ times
\begin{align*}
-\frac{\hbar^2}{2m}a''(x) - \frac{i\hbar}{m}a'(x)S'(x) - \frac{i\hbar}{2m}a(x)S''(x) + \frac{1}{2m}a(x)(S'(x))^2 + (V(x)-E)a(x).
\end{align*}
In the formal WKB sense, the leading equations are obtained by setting the coefficients of $\hbar^0$ and $\hbar^1$ in this expansion equal to zero; the remaining $\hbar^2$ term is a higher-order residual for the leading ansatz, not an equation imposed at this order. Thus the $\hbar^0$ coefficient is
\begin{align*}
a(x)\left(\frac{1}{2m}(S'(x))^2 + V(x) - E\right),
\end{align*}
and the $\hbar^1$ coefficient is
\begin{align*}
-\frac{i}{2m}\left(2a'(x)S'(x) + a(x)S''(x)\right).
\end{align*}
On a subinterval where $a$ is nonzero, the vanishing of the $\hbar^0$ coefficient permits division by $a(x)$ and gives
\begin{align*}
\frac{1}{2m}(S'(x))^2 + V(x) = E.
\end{align*}
The vanishing of the $\hbar^1$ coefficient gives
\begin{align*}
2a'(x)S'(x) + a(x)S''(x) = 0.
\end{align*}
[guided]
For each $x \in I$, define the Schrödinger residual map $R_\hbar:I\to\mathbb{C}$ by
\begin{align*}
R_\hbar(x) = -\frac{\hbar^2}{2m}\psi_\hbar''(x) + V(x)\psi_\hbar(x) - E\psi_\hbar(x).
\end{align*}
Substituting the second-derivative formula from the previous step and using $\psi_\hbar(x)=a(x)q_\hbar(x)$ gives
\begin{align*}
R_\hbar(x) = q_\hbar(x)\left[-\frac{\hbar^2}{2m}a''(x) - \frac{i\hbar}{m}a'(x)S'(x) - \frac{i\hbar}{2m}a(x)S''(x) + \frac{1}{2m}a(x)(S'(x))^2 + (V(x)-E)a(x)\right].
\end{align*}
The exponential factor satisfies $q_\hbar(x)=e^{iS(x)/\hbar}\neq 0$, so it cannot force the residual to vanish. The information is therefore contained in the bracketed formal expansion.
The phrase "coefficient gives" is interpreted in the formal WKB sense: for the leading ansatz, we set the coefficients of $\hbar^0$ and $\hbar^1$ equal to zero and regard the $\hbar^2$ term $-a''(x)/(2m)$ as belonging to the next order of the approximation. If instead one required this two-term expression to solve the Schrödinger equation exactly for every $\hbar>0$, then the $\hbar^2$ coefficient would also have to vanish; that is not the leading WKB convention used here.
The coefficient of $\hbar^0$ is
\begin{align*}
a(x)\left(\frac{1}{2m}(S'(x))^2 + V(x) - E\right).
\end{align*}
Setting this coefficient equal to zero gives a product equal to zero. Therefore, on any subinterval where the amplitude $a$ is nonzero, we may divide by $a(x)$ and obtain
\begin{align*}
\frac{1}{2m}(S'(x))^2 + V(x) = E.
\end{align*}
This is the Hamilton-Jacobi equation.
The coefficient of $\hbar^1$ is
\begin{align*}
-\frac{i}{2m}\left(2a'(x)S'(x) + a(x)S''(x)\right).
\end{align*}
Because $m>0$, the scalar $-i/(2m)$ is nonzero. Setting the $\hbar^1$ coefficient equal to zero is therefore equivalent to
\begin{align*}
2a'(x)S'(x) + a(x)S''(x) = 0.
\end{align*}
This is the transport equation for the leading amplitude.
[/guided]
[/step]
[step:Rewrite the transport equation as conservation of $a^2S'$]
Assume the transport equation holds on a subinterval $J \subset I$. No nonvanishing assumption on $S'$ is needed for the following conservation identity. Define $F: J \to \mathbb{C}$ by
\begin{align*}
F(x) = a(x)^2S'(x).
\end{align*}
Since $a \in C^2(I;\mathbb{C})$ and $S \in C^2(I;\mathbb{R})$, the function $F$ is $C^1$ on $J$. By the product rule,
\begin{align*}
F'(x) = 2a(x)a'(x)S'(x) + a(x)^2S''(x).
\end{align*}
Factoring out $a(x)$ gives
\begin{align*}
F'(x) = a(x)\left(2a'(x)S'(x) + a(x)S''(x)\right).
\end{align*}
The factor in parentheses is zero by the transport equation, so
\begin{align*}
(a(x)^2S'(x))' = 0.
\end{align*}
On any subinterval where $S'$ has no zero and $a$ is chosen nonzero, this conservation law is equivalently the usual leading amplitude relation $a(x)$ proportional to $|S'(x)|^{-1/2}$, up to a constant choice of branch.
[/step]
