Help deriving Lax-Wendroff scheme for advection equation $u_t+c(x)u_x = 0$

Question 1: Consider the wave equation $$u_t + c(x) u_x = 0 ,$$ where $x\in \Omega \subset \Bbb R$ and $c(x)$ is a function of $x$.

(a) Show that the Lax-Wendroff scheme for this PDE is given by $$u_j^{n+1} = u_j^n - c_j \Delta t \frac{D_x u_j^n}{2 \Delta x} + \frac{c_j^2 \Delta t^2}{2} \frac{\delta_x^2 u_j^n}{\Delta x^2} + \frac{c_j \Delta t^2}{8 \Delta x ^2} (D_x c_j)(D_x u_j^n) ,$$ where $D_x$ is the first central difference operator, $\delta_x^2$ is the second central difference operator, and $\Delta t$ and $\Delta x$ are the mesh-spacing in $t$ and $x$, respectively. The $j$ and $n$ are space and time indices, respectively, and $u_j^n$ is the grid function such that $u_j^n\approx u(x_j,t_n)$ and $c_j \approx c(x_j)$.

I need help deriving this particular scheme,

start by taylor expansion,

$u(x,t+\Delta t)= u(x,t) + u_t \Delta t +\ \frac{\Delta t^{2}}{2} u_{tt} + O (\Delta t)^{3}$

$u_t = \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t} -\frac{\Delta t^{2}}{2} u_{tt}$

$u_t = -cu_x$

$u_{tt} = c^{2} u_{xx}$

subing these in, $u_t = \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t} -\frac{\Delta t^{2}}{2} u_{xx}$

$u_t = \frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t} -\frac{ c^2 \Delta t^{2}}{ 2\Delta x^{2}} \delta x^{2} u_{j}^{n}$

$u _{x} =\frac{ u_{j+1}^{n}-u_{j-1}^{n}}{ 2 \Delta x }$ sub this in into $u_t+c(x)u_x = 0$

$u_{j}^{n+1} = u_{j}^{n} - \frac{1}{2}p u_{j+1}^{n}-u_{j-1}^{n} +\frac{1}{2} p^2( u_{j+1}^n-2u_{j}^{n} + u_{j-1}^{n} )$

Am really unsure how to derive this last part to the scheme, it seems so random, like its just put on or something?? Does anyone know how to derive this bit ?

We follow the steps in this post. The first and second-order time derivatives write $u_t = -c(x) u_x$ and \begin{aligned} u_{tt} &= -c(x) u_{tx} \\ & = c(x)\, (c(x) u_x)_x \\ &= c(x)\,( c'(x) u_x + c(x) u_{xx}) \, , \end{aligned} respectively. Using central finite differences in space, we therefore have $$u_t(x_j,t_n) \approx -c_j \frac{D_x u_j^n}{2 \Delta x}$$ and $$u_{tt}(x_j,t_n) \approx c_j \left(\frac{D_x c_j}{2 \Delta x}\frac{D_x u_j^n}{2\Delta x} + c_j \frac{\delta^2_x u_j^n}{\Delta x^2}\right) .$$ Injecting this Ansatz in the Taylor series $$u(x_j,t_{n+1}) = u(x_j,t_{n}) + \Delta t\, u_t(x_j,t_{n}) + \frac{\Delta t^2}{2} u_{tt}(x_j,t_{n}) +\dots$$ gives the proposed numerical method.