Similarity Transforms:

$\partial_t U_i(x,t) + \partial_{x_j} F_{ij}(U(x,t)) = 0$

substitute $\tilde{F}_{ij}$ and $\tilde{U}_i$ into conservative equations
$\partial_t \tilde{U}_i + \partial_{x_j} (\tilde{F}_{ij} + \xi_j \tilde{U}_i) = 0$

time derivative of state (in self-similar form)
$\partial_t \tilde{U}_i = {\partial\tilde{U}_i\over\partial\xi_j} {\partial\xi_j\over\partial t} = -{x_j\over t^2} \partial_{\xi_j} \tilde{U}_i$
substitute $\xi_j = {x_j\over t}$
$\partial_t \tilde{U}_i = -{1\over t} \xi_j \partial_{\xi_j} \tilde{U}_i$

spatial derivative of flux (in self-similar form)
$\partial_{x_j} F_{ij} = \partial_{x_j} (\tilde{F}_{ij} + \xi_j \tilde{U}_i)$
$= {{\partial \xi_j}\over{\partial x_j}} \partial_{\xi_j} (\tilde{F}_{ij} + \xi_j \tilde{U}_i)$
$= {1\over t} (\partial_{\xi_j} \tilde{F}_{ij} + (\partial_{\xi_j} \xi_j) \tilde{U}_i + \xi_j \partial_{\xi_j} \tilde{U}_i)$
$= {1\over t} (\partial_{\xi_j} \tilde{F}_{ij} + n \tilde{U}_i) + {1\over t} \xi_j \partial_{\xi_j} \tilde{U}_i$ for $n$ dimensions
substitute to find:
$ -{1\over t} \xi_j \partial_{\xi_j} \tilde{U}_i + {1\over t} (\partial_{\xi_j} \tilde{F}_{ij} + n \tilde{U}_i) + {1\over t} \xi_j \partial_{\xi_j} \tilde{U}_i = 0$

simplifies to
$n \tilde{U}_i + \partial_{\xi_j} \tilde{F}_{ij} = 0$

finite volume form:
$\int_\Omega (n \tilde{U}_i + \partial_{\xi_j} \tilde{F}_{ij}) dx = 0$
$\int_\Omega dx n \tilde{U}_i + \int_\Omega \partial_{\xi_j} \tilde{F}_{ij} dx = 0$
$|\Omega| n \tilde{U}_i + \int_{\partial\Omega} \tilde{F}_{ij} n_j d\sigma = 0$
discretize:
$\Delta \xi n \tilde{U}_i + \tilde{F}_{ij}(x_R) - \tilde{F}_{ij}(x_L) = 0$
$n \tilde{U}_i + {1\over\Delta \xi}(\tilde{F}_{ij}(x_R) - \tilde{F}_{ij}(x_L)) = 0$
...for $\Delta \xi = \Pi_{k=1}^n \Delta \xi_k$

Similarity Transform of MHD Equations:

$U_i = \tilde{U}_i = \left[\begin{matrix} \rho \\ \rho v_i \\ B_i \\ E_{total} \end{matrix}\right]$

$F_{ij} = \left[\begin{matrix} \rho v_j \\ \rho v_i v_j + g^{ij} P_{total} - \frac{1}{\mu_0} B_i B_j \\ B_i v_j - B_j v_i \\ (E_{total} + P_{total}) v_j - \frac{1}{\mu_0} B_k v_k B_j \end{matrix}\right]$

$\tilde{F}_{ij} = F_{ij} - \xi_j \tilde{U}_i = \left[\begin{matrix} \rho (v_j - \xi_j) \\ \rho v_i (v_j - \xi_j) + g^{ij} P_{total} - \frac{1}{\mu_0} B_i B_j \\ B_i (v_j - \xi_j) - B_j v_i \\ E_{total} (v_j - \xi_j) + P_{total} v_j - \frac{1}{\mu_0} B_k v_k B_j \end{matrix}\right]$

Newton Method

Nonlinear function:
$G_i(\tilde{U}) = n \tilde{U}_i + \partial_{\xi_j} \tilde{F}_{ij}$
solve for $G_i = 0$ by minimizing $\tilde{U}_j$

Taylor expansion of function to minimize:
$G_i(\tilde{U}^{k+1}) = G_i(\tilde{U}^k) + {{\partial G_i(\tilde{U}^k)}\over{\partial \tilde{U}^k_j}} (\tilde{U}_j^{k+1} - \tilde{U}_j^k) + \mathscr(O)((\tilde{U})^2)$
$\tilde{U}_i^k$ the $k$th iteration of our state vector $\tilde{U}_i$
$G_i(\tilde{U}^{k+1}) = G_i(\tilde{U}^k) + {{\partial G_i(\tilde{U}^k)}\over{\partial \tilde{U}^k_j}} (\tilde{U}_j^{k+1} - \tilde{U}_j^k) + \mathscr(O)((\tilde{U})^2)$

Rearranged, and ignoring the higher-order terms:
$||{{\partial G_i(\tilde{U}^k)}\over{\partial \tilde{U}^k_j}}||^{-1}_{ij} (G_j(\tilde{U}^{k+1}) - G_j(\tilde{U}^k)) = \tilde{U}_i^{k+1} - \tilde{U}_i^k$
$\tilde{U}_i^{k+1} = \tilde{U}_i^k + ||{{\partial G_i(\tilde{U}^k)}\over{\partial \tilde{U}^k_j}}||^{-1}_{ij} (G_j(\tilde{U}^{k+1}) - G_j(\tilde{U}^k))$
Assume that we arrive at our destination, i.e. $G_i(\tilde{U}^{k+1}) = 0$:
$\tilde{U}_i^{k+1} = \tilde{U}_i^k - ||{{\partial G_i(\tilde{U}^k)}\over{\partial \tilde{U}^k_j}}||^{-1}_{ij} G_j(\tilde{U}^k)$
This formula is equivalent to Newton's method of an n-dimensions scalar function with value equal to $\frac{1}{2}||G||^2$

Jacobian Approximation:
${{\partial G_i(u)}\over{\partial \tilde{U}_j}} v_j \approx ( G_i (u + \epsilon v) - G_i(u) ) / \epsilon$

Evaluation of function to minimize
${{\partial G_i}\over{\partial \tilde{U}_j}} = \partial_{\tilde{U}_j} G_i$
$= \partial_{\tilde{U}_j} (n \tilde{U}_i + \partial_{\xi_k} \tilde{F}_{ik})$
$= n g^{ij} + \partial_{\tilde{U}_j} \partial_{\xi_k} \tilde{F}_{ik}$