Lax-Wendroff method for nonlinear hyperbolic systems of conservation laws In page 127 of R.J. LeVeque's "Numerical Methods for conservation laws" (Birkhäuser, 1992),

There are various ways that [the Lax-Wendroff method for constant-coefficient linear hyperbolic systems] can be extended to give a second order method for nonlinear conservation laws. If we let $A(u) = f'(u)$ be the Jacobian matrix, then a conservative form of Lax-Wendroff is
  $$
U_j^{n+1} = U_j^n - \frac{k}{2h}(f(U_{j+1}^n - f(U_{j-1}^n)) + \frac{k^2}{2h^2} \big[ A_{j+1/2} (f(U_{j+1}^n) - f(U_j^n)) \\
- A_{j-1/2} (f(U_{j}^n) - f(U_{j-1}^n)) \big] \, ,
\tag{12.25}
$$
  where $A_{j\pm1/2}$ is the Jacobian matrix evaluated at $\frac12 (U_j^n + U_{j\pm 1}^n)$.

I am trying to understand the paragraph above, which is pretty much the same as in the Wikipedia article. First of all, how was this derived and how is it in conservative form? Moreover, Why are we evaluating a jacobian at $\frac12 (U_j^n + U_{j \pm 1}^n )$? Isnt this the value of the solution $u= u(x)$ evaluated at some discrete point? Im confused by this notation. Can someone clarify this? 
My question is related to an exercise of this book:

EXERCISE 12.3. [The method] $(12.25)$ reduces to [the classical linear Lax-Wendroff method] in the constant coefficient linear case and is second order accurate on smooth solutions (to nonlinear problems) and conservative. Verify these statements for at least [$(12.25)$] and write it in conservation form, determining the numerical flux function.

 A: Start with the Taylor expansion 
$$
u(x,t+k)=u(x,t)+u_t(x,t)k+\frac12u_{tt}(x,t)k^2+...
$$
Use the PDE to convert time into space derivatives, $u_t=-[f(u)]_x$, then
\begin{align}
u_{tt}&=-[f(u)]_{xt}=-[f'(u)u_t]_x\\
&=[f'(u)[f(u)]_x]_x\\
\end{align}
Now realize this derivative structure using symmetric divided differences and a half-step scheme. Let the index denote the offset in $x$ direction in units of $h$
$$
[f'(u)[f(u)]_x]_x=\frac{f'(u_{+1/2})[f(u_{+1/2})]_x-f'(u_{-1/2})[f(u_{-1/2})]_x}h
$$
The $x$ derivatives at the half-steps are again approximated using central divided differences with step size $h/2$, 
$$
[f'(u)[f(u)]_x]_x=\frac{f'(u_{+1/2})\frac{f(u_{+1})-f(u)}h-f'(u_{-1/2})\frac{f(u)-f(u_{-1})}h}h
$$
and finally the midpoint arguments of the derivative of $f$ by $u$ can be approximated by the mean of the neighboring integer grid values, $u_{\pm 1/2}=\frac12(u+u_{\pm 1})$.
If one inserts this in backwards, using the full-step central derivative for the degree-1 term, $[f(u)]_x=\frac{f(u_{+1})-f(u_{-1})}{2h}$ one gets back to the given formula of the method.
A: Just one more comment, which you probably know. The problem of the Lax-Wendroff scheme is that it contains zero numerical viscosity and therefore it cannot filter out oscillation which the scheme produces due to its second order approximation. For the practical use I recommend the TVD or ENO schemes, which are more sofisticated.
