Editing Appendix/Ramblings/51BiPolytropeStability/BetterInterface (section)

==Discretize for Numerical Integration==

===General Discretization===

====First Approximation====
Now, let's set up a grid associated with a uniformly spaced spherical radius, where the subscript <math>J</math> denotes the grid zone at which all terms in the finite-difference representation of the governing relations will be evaluated.  More specifically,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\tilde{r}_{J-1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\tilde{r}_J - \Delta\tilde{r}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\tilde{r}_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\tilde{r}_J + \Delta\tilde{r} \, ;
</math>
  </td>
</tr>
</table>
also,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\biggl(\frac{dx}{d\tilde{r}}\biggr)_{J}</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
\frac{(x_{J+1} - x_{J-1})}{2\Delta\tilde{r}}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\biggl(\frac{dp}{d\tilde{r}}\biggr)_{J}</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
\frac{(p_{J+1} - p_{J-1})}{2\Delta\tilde{r}} \, .
</math>
  </td>
</tr>
</table>
And at each grid location, the governing relations establish the local evaluation of the derivatives, that is,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\biggl(\frac{dx}{d \tilde{r}}\biggr)_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
- \frac{1}{\tilde{r}_J}\biggl[ 
3x + \frac{p}{\gamma_g}\biggr]_J \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\biggl(\frac{dp}{d \tilde{r}}\biggr)_J </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{\tilde{\rho}_J}{\tilde{P}_J}\biggl[ 
(4x + p)\frac{\tilde{M}_r}{\tilde{r}^2}
+
\tau_c^2 \omega^2 \tilde{r} x \biggr]_J \, .
</math>
  </td>
</tr>
</table>
<span id="1stapprox">So, integrating</span> step-by-step from the center of the configuration, outward, once all the variable values are known at grid locations <math>J</math> and <math>(J-1)</math>, the values of <math>x</math> and <math>p</math> at <math>(J+1)</math> are given by the expressions,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} - 2\Delta\tilde{r} \biggl\{
\frac{1}{\tilde{r}}\biggl[ 
3x + \frac{p}{\gamma_g}\biggr]
\biggr\}_J \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>p_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + 2\Delta\tilde{r}
\biggl\{
\frac{\tilde{\rho}}{\tilde{P}} \cdot \frac{\tilde{M}_r}{\tilde{r}^2}
\biggl[ (4x + p)
+
\sigma_c^2 \biggl(\frac{2\pi}{3} \cdot \frac{\tilde{\rho}_c \tilde{r}^3 }{\tilde{M}_r}\biggr) x \biggr] 
\biggr\}_J\, .
</math>
  </td>
</tr>
</table>

Then we will obtain the "<math>x_J</math>" and "<math>p_J</math>" values via the ''average'' expressions,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{2}(x_{J-1} + x_{J+1})
\, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>p_{J}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{2}(p_{J-1} + p_{J+1})
\, .
</math>
  </td>
</tr>
</table>

<table border="1" align="center" cellpadding="8" width="80%"><tr><td align="left">
Consider implementing a more ''implicit'' finite-difference analysis.  Wherever a "<math>J</math>" index appears in the source term, replace it with the ''average expressions.''  The general form of the source term expressions is,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + 2\Delta\tilde{r} \biggl\{
\mathcal{A}x_J + \mathcal{B}p_J
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + \Delta\tilde{r} \biggl\{
\mathcal{A}\biggl[x_{J+1}+x_{J-1}\biggr] + \mathcal{B}\biggl[p_{J+1}+p_{J-1}\biggr]
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ x_{J+1}\biggl[1 - \Delta\tilde{r} \mathcal{A}  \biggr]</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{A}\biggr] + \Delta\tilde{r} 
\mathcal{B}\biggl[p_{J+1}+p_{J-1}\biggr]
\, ,
</math>
  </td>
</tr>
</table>
where,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\mathcal{A}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
- \frac{3}{\tilde{r}} \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\mathcal{B}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
- \frac{1}{\gamma_g \tilde{r}} \, ;
</math>
  </td>
</tr>
</table>
and,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + 2\Delta\tilde{r} \biggl\{
\mathcal{C}x_J + \mathcal{D}p_J
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1} + \Delta\tilde{r} \biggl\{
\mathcal{C}\biggl[x_{J+1}+x_{J-1}\biggr] + \mathcal{D}\biggl[p_{J+1}+p_{J-1}\biggr]
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ p_{J+1}\biggl[1 - \Delta\tilde{r} \mathcal{D}  \biggr]</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
p_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{D}\biggr] + \Delta\tilde{r} 
\mathcal{C}\biggl[x_{J+1}+x_{J-1}\biggr]
\, ,
</math>
  </td>
</tr>
</table>

where,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>\mathcal{C}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
\biggl\{
\mathcal{D}
\biggl[ 4
+
\sigma_c^2 \biggl(\frac{2\pi}{3} \cdot \frac{\tilde{\rho}_c \tilde{r}^3 }{\tilde{M}_r}\biggr)  \biggr] 
\biggr\}_J\, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>\mathcal{D}</math></td>
  <td align="center"><math>\equiv</math></td>
  <td align="left">
<math>
\biggl\{
\frac{\tilde{\rho}}{\tilde{P}} \cdot \frac{\tilde{M}_r}{\tilde{r}^2}
\biggr\}_J\, .
</math>
  </td>
</tr>
</table>

In both cases, the two unknowns are <math>x_{J+1}</math> and <math>p_{J+1}</math>.  Combining this pair of equations gives,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}\biggl[1 - \Delta\tilde{r} \mathcal{A}  \biggr]</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{A}\biggr] 
+ \Delta\tilde{r} \mathcal{B}\biggl[p_{J-1}\biggr]
+ \Delta\tilde{r} \mathcal{B}\biggl[1 - \Delta\tilde{r} \mathcal{D}  \biggr]^{-1}\biggl\{ 
p_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{D}\biggr] + \Delta\tilde{r} 
\mathcal{C}\biggl[x_{J+1}+x_{J-1}\biggr]
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ x_{J+1}\biggl\{ 1 - \Delta\tilde{r} \mathcal{A}  
- \Delta\tilde{r} \mathcal{B}\biggl[1 - \Delta\tilde{r} \mathcal{D}  \biggr]^{-1}\biggl[ \Delta\tilde{r} \mathcal{C}\biggr]
\biggr\}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{A}\biggr] 
+ \Delta\tilde{r} \mathcal{B}\biggl[p_{J-1}\biggr]
+ \Delta\tilde{r} \mathcal{B}\biggl[1 - \Delta\tilde{r} \mathcal{D}  \biggr]^{-1}\biggl\{ 
p_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{D}\biggr] + \Delta\tilde{r} 
\mathcal{C}\biggl[x_{J-1}\biggr]
\biggr\} \, ,
</math>
  </td>
</tr>
</table>
which determines <math>x_{J+1}</math>, which then allows the straightforward determination of <math>p_{J+1}</math>.  Via the ''average'' expressions, we can also then determine &#8212; and record &#8212; the self-consistent values of <math>x_J</math> and <math>p_J</math>.

----

Dropping terms <math>\mathcal{O}[(\Delta\tilde{r})^2]</math> and higher gives,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}\biggl[ 1 - \Delta\tilde{r} \mathcal{A}  
\biggr]_J</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
x_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{A}\biggr]_J 
+ \Delta\tilde{r} \mathcal{B}_J\biggl[p_{J-1}\biggr]
\, ,
</math>
  </td>
</tr>
</table>
and, in turn,
<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>p_{J+1}\biggl[1 - \Delta\tilde{r} \mathcal{D}  \biggr]_J</math></td>
  <td align="center"><math>\approx</math></td>
  <td align="left">
<math>
p_{J-1}\biggl[1 + \Delta\tilde{r}\mathcal{D}\biggr]_J 
+ \Delta\tilde{r} \mathcal{C}_J\biggl[x_{J-1}\biggr]
+ \Delta\tilde{r} \mathcal{C}_J\biggl[
x_{J+1}
\biggr] \, .
</math>
  </td>
</tr>
</table>

</td></tr></table>

====Second Approximation====

Let's assume that we know the three quantities, <math>x_{J-1}, x_J</math>, and <math>(x_J)^' \equiv (dx/d\tilde{r})_J</math> and want to project forward to determine, <math>x_{J+1}</math>.  We should assume that, locally, the displacement function <math>x</math> is quadratic in <math>\tilde{r}</math>, that is,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b\tilde{r} + c\tilde{r}^2 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ \frac{dx}{d\tilde{r}}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
b + 2c\tilde{r} \, ,
</math>
  </td>
</tr>
</table>
where we have three unknowns, <math>a, b, c</math>.  These can be determined by appropriately combining the three relations,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>(x_J)^'</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
b + 2c\tilde{r}_J \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>x_J</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b\tilde{r}_J + c\tilde{r}_J^2 \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>x_{J-1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b(\tilde{r}_{J}-\Delta\tilde{r}) + c(\tilde{r}_{J}-\Delta\tilde{r})^2 \, .
</math>
  </td>
</tr>
</table>
We have,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>b</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
(x_J)^' - 2c\tilde{r}_J \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>- a</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-x_J + [(x_J)^' - 2c\tilde{r}_J]\tilde{r}_J + c\tilde{r}_J^2 \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>- a</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-x_{J-1} + [(x_J)^' - 2c\tilde{r}_J](\tilde{r}_{J}-\Delta\tilde{r}) + c(\tilde{r}_{J}-\Delta\tilde{r})^2 \, .
</math>
  </td>
</tr>
</table>
Combining the last two expressions gives,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>-x_J + [(x_J)^' - 2c\tilde{r}_J]\tilde{r}_J + c\tilde{r}_J^2</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-x_{J-1} + [(x_J)^' - 2c\tilde{r}_J](\tilde{r}_{J}-\Delta\tilde{r}) + c(\tilde{r}_{J}-\Delta\tilde{r})^2
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ -x_J + (x_J)^'\tilde{r}_J - c\tilde{r}_J^2 </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-x_{J-1} 
+ (x_J)^'(\tilde{r}_{J}-\Delta\tilde{r}) 
- 2c\tilde{r}_J(\tilde{r}_{J}-\Delta\tilde{r}) 
+ c[ \tilde{r}_{J}^2 - 2r_J \Delta\tilde{r} + (\Delta\tilde{r})^2 ]
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ -x_J  - c\tilde{r}_J^2 </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-x_{J-1} 
- (x_J)^'\Delta\tilde{r} 
+ c[ \tilde{r}_{J}^2 - 2r_J \Delta\tilde{r} + (\Delta\tilde{r})^2 
- 2\tilde{r}_J^2 + 2\tilde{r}_J(\Delta\tilde{r})]
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ c(\Delta\tilde{r})^2
</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J -x_{J-1} - (x_J)^'\Delta\tilde{r}
\, .
</math>
  </td>
</tr>
</table>
Therefore, also,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>b</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
(x_J)^' - \biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{2\tilde{r}_J}{(\Delta\tilde{r})^2}
 \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>a</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^' \tilde{r}_J 
+ \biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{\tilde{r}_J^2}{(\Delta\tilde{r})^2}
 </math>
  </td>
</tr>
</table>

Hence,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b(\tilde{r}_J + \Delta\tilde{r}) + c(\tilde{r}_J + \Delta\tilde{r})^2 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^' \tilde{r}_J 
+ \biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{\tilde{r}_J^2}{(\Delta\tilde{r})^2}
+ 
(\tilde{r}_J + \Delta\tilde{r})\biggl\{
(x_J)^' - \biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{2\tilde{r}_J}{(\Delta\tilde{r})^2}
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center">&nbsp;</td>
  <td align="left">
<math>
+ 
(\tilde{r}_J + \Delta\tilde{r})^2 
\biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{1}{(\Delta\tilde{r})^2}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^' \tilde{r}_J 
+ 
\biggl\{
\tilde{r}_J^2 
+ 2\tilde{r}_J(\tilde{r}_J + \Delta\tilde{r})
+ (\tilde{r}_J + \Delta\tilde{r})^2
\biggr\} 
\biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{1}{(\Delta\tilde{r})^2}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^' \tilde{r}_J 
+ 
\biggl[ 2\tilde{r}_J + \Delta\tilde{r}\biggr]^2
\biggl[ x_J -x_{J-1} - (x_J)^'\Delta\tilde{r} \biggr]\frac{1}{(\Delta\tilde{r})^2}
\, .
</math>
  </td>
</tr>
</table>
<font color="red"><b>WRONG!!</b></font>  Try again &hellip;

====Third Approximation====

Let's assume that we know the three quantities, <math>x_{J-1}, x_J</math>, and <math>(x_J)^' \equiv (dx/d\tilde{r})_J</math> and want to project forward to determine, <math>x_{J+1}</math>.  We should assume that, locally, the displacement function <math>x</math> is quadratic in <math>\tilde{r}</math>, that is,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b\tilde{r} + c\tilde{r}^2 
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\Rightarrow ~~~ \frac{dx}{d\tilde{r}}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
b + 2c\tilde{r} \, ,
</math>
  </td>
</tr>
</table>
where we have three unknowns, <math>a, b, c</math>.  These can be determined by appropriately combining the three relations,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>(x_J)^'</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
b + 2c\tilde{r}_J \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>x_J</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b\tilde{r}_J + c\tilde{r}_J^2 \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>x_{J-1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
a + b(\tilde{r}_{J}-\Delta\tilde{r}) + c(\tilde{r}_{J}-\Delta\tilde{r})^2 \, .
</math>
  </td>
</tr>
</table>
The difference between the last two expressions gives,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_J - x_{J-1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl[ b\tilde{r}_J + c\tilde{r}_J^2 \biggr]
-
\biggl[ b(\tilde{r}_{J}-\Delta\tilde{r}) + c(\tilde{r}_{J}-\Delta\tilde{r})^2 \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
b\Delta\tilde{r} 
+ 
c(2\tilde{r}_J \Delta\tilde{r} - \Delta\tilde{r}^2) \, .
</math>
  </td>
</tr>
</table>
Combining this with the first of the three expressions then gives,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right">
<math>
\Delta\tilde{r} \biggl[(x_J)^' - 2c\tilde{r}_J  \biggr]
+ 
c(2\tilde{r}_J \Delta\tilde{r} - \Delta\tilde{r}^2) 
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - x_{J-1}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ 
c(2\tilde{r}_J \Delta\tilde{r} - \Delta\tilde{r}^2) 
-c\biggl[ 2\tilde{r}_J \Delta\tilde{r}  \biggr]
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ 
-c \Delta\tilde{r}^2 
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ 
c  
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{\Delta\tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\, .
</math>
  </td>
</tr>
</table>
Hence,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right">
<math>
b
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
(x_J)^' - 2c\tilde{r}_J 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
(x_J)^' + \frac{2\tilde{r}_J}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] 
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
a
</math>
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - b\tilde{r}_J -c\tilde{r}_J^2
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J
-\biggl\{ (x_J)^' + \frac{2\tilde{r}_J}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] \biggr\}\tilde{r}_J
- \biggl\{
\frac{1}{\Delta\tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}\tilde{r}_J^2
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^'\tilde{r}_J
-\biggl\{  \frac{2\tilde{r}_J^2}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] \biggr\}
+ \biggl\{
\frac{\tilde{r}_J^2}{\Delta\tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^'\tilde{r}_J
-\frac{\tilde{r}_J^2}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] 
\, .
</math>
  </td>
</tr>
</table>
As a result,

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{
a 
\biggr\}
+ (\tilde{r}_J +\Delta\tilde{r}) \biggl\{
b
\biggr\}
+ (\tilde{r}_J+\Delta\tilde{r})^2 \biggl\{
c
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl\{
x_J - (x_J)^'\tilde{r}_J
-\frac{\tilde{r}_J^2}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
+ (\tilde{r}_J +\Delta\tilde{r}) \biggl\{
(x_J)^' + \frac{2\tilde{r}_J}{\Delta \tilde{r}^2}\biggl[
x_J - x_{J-1} - (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center">&nbsp;</td>
  <td align="left">
<math>
+ (\tilde{r}_J^2 + 2\tilde{r}_J \Delta\tilde{r} +\Delta\tilde{r}^2) \biggl\{
\frac{1}{\Delta\tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J - (x_J)^'\tilde{r}_J
+ \tilde{r}_J^2\biggl\{\frac{1}{\Delta \tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
+ (x_J)^'(\tilde{r}_J +\Delta\tilde{r}) 

- 2\tilde{r}_J(\tilde{r}_J +\Delta\tilde{r}) \biggl\{
\frac{1}{\Delta \tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center">&nbsp;</td>
  <td align="left">
<math>
+ (\tilde{r}_J^2 + 2\tilde{r}_J \Delta\tilde{r} +\Delta\tilde{r}^2) \biggl\{
\frac{1}{\Delta\tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J + (x_J)^'\Delta\tilde{r} 

+ \biggl[ \tilde{r}_J^2 - 2\tilde{r}_J(\tilde{r}_J +\Delta\tilde{r})+ (\tilde{r}_J^2 + 2\tilde{r}_J \Delta\tilde{r} +\Delta\tilde{r}^2)\biggr]
\biggl\{\frac{1}{\Delta \tilde{r}^2}\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_J + (x_J)^'\Delta\tilde{r} 
+ 
\biggl[
-x_J + x_{J-1} + (x_J)^'\Delta\tilde{r}
\biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + 2(x_J)^'\Delta\tilde{r} \, .
</math>
  </td>
</tr>
</table>
<font color="red"><b>GOOD!</b></font>  This is the same as our [[#1stapprox|first approximation expression]] stated above.


<table border=1 align="center" cellpadding="10" width="80%"><tr><td bgcolor="lightgreen" align="left">
This is test ...
<table border="1" align="center" cellpadding="5">

  <tr>
<td align="center" bgcolor="white"><math>\Delta\tilde{r}</math></td>
<td align="center" bgcolor="white"><math>x_J</math></td>
<td align="center" bgcolor="white"><math>x_{J-1}</math></td>
<td align="center" bgcolor="white"><math>(x_J)^'</math></td>
<td align="center" bgcolor="white"><math>(x_{J-1})^'</math></td>
  </tr>

  <tr>
<td align="center" bgcolor="white">0.001936393</td>
<td align="center" bgcolor="white">-4.695376</td>
<td align="center" bgcolor="white">-4.547832</td>
<td align="center" bgcolor="white">-116.0119</td>
<td align="center" bgcolor="white">-76.19513</td>
  </tr>
</table>
<tr><td bgcolor="white" align="center">

<table border="0" align="center" cellpadding="5">

<tr>
  <td align="right"><math>x_{J+1}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
x_{J-1} + 2(x_J)^'\Delta\tilde{r}
=
-4.997121
\, .
</math>
  </td>
</tr>
</table>
</td></tr>
</td></tr></table>