Editing SSC/Structure/BiPolytropes/AnalyzeStepFunction (section)

=Discrete LAWE for Bipolytrope=

==Summary Set of Nonlinear Governing Relations==
In summary, the following three one-dimensional ODEs define the physical relationship between the three dependent variables {{Math/VAR_Density01}}, {{Math/VAR_Pressure01}}, and <math>~r</math>, each of which should be expressible as a function of the two independent (Lagrangian) variables, {{Math/VAR_Time01}} and <math>~M_r</math>:

<div align="center">
<span id="Continuity"><font color="#770000">'''Equation of Continuity'''</font></span><br />
<math>\frac{d\rho}{dt} = - 4\pi \rho^2 r^2  \frac{d}{dM_r}\biggl(\frac{dr}{dt}\biggr) - \frac{2\rho}{r} \biggl(\frac{dr}{dt}\biggr)  </math><br /> ,

<span id="PGE:Euler"><font color="#770000">'''Euler + Poisson Equations'''</font></span><br />
<math>\frac{d^2 r}{dt^2} = - 4\pi r^2 \frac{dP}{dM_r} - \frac{GM_r}{r^2} </math><br />


<span id="PGE:AdiabaticFirstLaw">Adiabatic Form of the<br />
<font color="#770000">'''First Law of Thermodynamics'''</font></span><br />
<math>
\rho \frac{dP}{dt} - \gamma_\mathrm{g} P \frac{d\rho}{dt} = 0 .
</math>

</div>

==March from the Center, Outward==

===Establish Fidelity of Finite-Difference Model===
We know the analytic structure of the equilibrium configuration.  Let's choose a Lagrangian grid that is labeled by <math>(r_0)_j</math> and the corresponding enclosed mass, <math>m_j(r_0)</math>, where the center of the spherical bipolytrope is labeled by <math>j=0</math> while each subsequent grid "line" is labeled <math>j</math>.  We will identify the mean density of each mass shell by the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>{\bar\rho}_{j-1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[\frac{\Delta m}{\mathrm{Volume}}\biggr]_{j-1/2}
=
\biggl[m_j - m_{j-1}\biggr]\biggl[\frac{4\pi}{3} (r_j^3 - r_{j-1}^3)  \biggr]^{-1} \, .
</math>
  </td>
</tr>
</table>
The pressure can be determined from knowledge of the density via knowledge of the (fixed) specific entropy, namely,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>P_{j-1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
({\bar\rho}_{j-1/2})^{\gamma} \cdot \exp\biggl[ \frac{\mu (\gamma-1)s}{\mathfrak{R}} \biggr]
\, .
</math>
  </td>
</tr>
</table>
These two expressions, effectively, originate from the continuity equation and the adiabatic form of the first law of thermodynamics, respectively.  They are relations that allow the determination of the mass density and the pressure, given fixed mass shells but varying mass-shell radial locations.

<font color="red"><b>STEP 1:</b></font> &nbsp; From the analytically known equilibrium structure of the <math>(n_c, n_e) = (5, 1)</math> bipolytrope, create a table that documents how the radial location of each mass shell, <math>r_j</math>, varies with the enclosed mass, <math>m_j</math>.

<font color="red"><b>STEP 2:</b></font> &nbsp; Determine how <math>{\bar\rho}_{j-1/2}</math> varies with radial shell location, using the above continuity-equation relation.  (Plot <math>{\bar\rho}</math> versus <math>m</math> obtained in this discrete manner &#8212; see the red-dotted curve in the following figure &#8212; then also plot how <math>\rho</math> varies with <math>m</math> according to the analytic equilibrium structure &#8212; see the dark-blue-dotted curve; the plotted curves should be nearly, but not exactly, the same.) 

<font color="red"><b>STEP 3A:</b></font> &nbsp; Given <math>{\bar\rho}_{j-1/2}</math>, determine how <math>{P}_{j-1/2}</math> varies with radial shell location, using the above 1<sup>st</sup> Law relation.  (Plot the pressure, as determined in this discrete manner, versus <math>m</math> &#8212; see the light-green-dotted curve in the following figure &#8212; then also plot how <math>P</math> varies with <math>m</math> according to the analytic equilibrium structure &#8212; see the purple-dotted curve; the plotted curves should be nearly, but not exactly, the same.)

<font color="red"><b>STEP 3B:</b></font> &nbsp; Determine (and plot) how <math>(dP/dm)_{j}</math> varies with <math>m</math>.  

Now, what can we learn from the "Euler + Poisson Equation"?  Well, for the equilibrium state, we should find that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{dP}{dm}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{Gm}{4\pi r^4}
\, .
</math>
  </td>
</tr>
</table>

<font color="red"><b>STEP 4:</b></font> &nbsp; Show how <math>dP/dm</math> varies with <math>m</math>, according to this analytic prescription, and compare it against the pressure gradient behavior obtained in STEP 3.  Do they match?

<table border="1" width="80%" align="center" cellpadding="5">
<tr>
  <td align="center">[[File:FiniteDifferenceLAWE01.png|600px|Finite-Difference Representation of n = 5 Core]]</td>
</tr>
<tr>
  <td align="left">
On a semi-log plot &hellip;
* As prescribed in <font color="red">STEP 2</font>: The red-dotted curve shows how the ''discrete'' density, <math>{\bar\rho}_{j-1/2}</math>, varies with enclosed mass, and the dark-blue-dotted curve shows how the analytically determined density varies with enclosed mass. 
* As prescribed in <font color="red">STEP 3A</font>: The light-green-dotted curve shows how the ''discrete'' pressure, <math>{P}_{j-1/2}</math>, varies with enclosed mass, and the purple-dotted curve shows how the analytically determined pressure varies with enclosed mass. 
* As prescribed in <font color="red">STEP 3B</font>: The light-blue-dotted curve shows how the ''discrete'' pressure gradient, <math>(dP/dm)_j</math>, varies with enclosed mass, and the orange-dotted curve shows how the analytically determined quantity, <math>-Gm/(4\pi r^4)</math>, varies with enclosed mass. Note that, in each case, we have added "1" to the logarithm of the specified quantity in order to shift both curves up in the plot and thereby unclutter the diagram.
  </td>
</tr>
</table>

===Guessing Eigenvector===

====Step 5====
<font color="red"><b>STEP 5:</b></font> &nbsp; Guess the eigenvector, <math>{\delta r}_i</math>, remembering that a reasonably good trial eigenfunction for the core is one that has a "parabolic dependence on the radius," namely,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{x}{\alpha_\mathrm{scale}}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
1 - \frac{\xi^2}{15} \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; where, &nbsp; &nbsp; &nbsp; 
  <td align="right">
<math>\xi</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(\frac{2\pi}{3}\biggr)^{1 / 2} r_0 \, .
</math>
  </td>
</tr>
</table>

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

<tr>
  <td align="right">
<math>x(r_0)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale}\biggl[1 - \biggl(\frac{2\pi}{45}\biggr)r_0^2  \biggr]
\, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \frac{dx}{dr_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
-\alpha_\mathrm{scale}\biggl(\frac{4\pi}{45}\biggr)r_0 
\, .
</math>
  </td>
</tr>
</table>

Note as well that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{d\ln x}{d\ln r_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
-\alpha_\mathrm{scale}\biggl(\frac{4\pi}{45}\biggr) \frac{r_0^2}{x}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(\frac{4\pi}{45}\biggr) \frac{3\xi^2}{2\pi} \biggl[\frac{\xi^2}{15}-1\biggr]^{-1}
\, .
</math>
  </td>
</tr>
</table>

<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<div align="center"><b>Finite-Difference Representation of <math>\delta r</math></b></div>
Once a "guess" for the fractional displacement vector, <math>(\delta r)_j = [x \cdot r_0]_j\, ,</math> has been specified, we recognize that the perturbed location of each radial shell is given by the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>r_j</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
(r_0)_j + (\delta r)_j \, .
</math>
  </td>
</tr>
</table>
</td></tr></table>


----


<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center" colspan="13">Parabolic Displacement Function w/ <math>\alpha_\mathrm{scale} = -0.001</math></td>
</tr>

<tr>
  <td align="center" rowspan="2">Shell</td>
  <td align="center" rowspan="2"><math>m </math></td>
  <td align="center" rowspan="2"><math>r_0</math></td>
  <td align="center" rowspan="2"><math>\frac{x}{\alpha_\mathrm{scale}}</math></td>
  <td align="center" rowspan="2"><math>r_j</math></td>
  <td align="center" rowspan="1" colspan="2">Analytic</td>
  <td align="center" rowspan="1" colspan="3">FD</td>
  <td align="center" rowspan="2"><math>P^*</math></td>
  <td align="center" rowspan="2"><math>-~\frac{dP^*}{dm}</math></td>
  <td align="center" rowspan="2"><math>-4\pi (r^*)^2\frac{dP^*}{dm} = g_0</math></td>
</tr>

<tr>
  <td align="center" rowspan="1"><math>\rho_0</math></td>
  <td align="center" rowspan="1"><math>d_\mathrm{analytic}</math></td>
  <td align="center" rowspan="1"><math>[\rho_{j - 1 / 2}]_0</math></td>
  <td align="center" rowspan="1"><math>d_{j-1 / 2}</math></td>
  <td align="center" rowspan="1"><math>\rho_{j - 1 / 2}</math></td>
</tr>

<tr>
  <td align="center">0</td>
  <td align="right">0.000000</td>
  <td align="right">0.000000</td>
  <td align="right">1.000000</td>
  <td align="right">0.000000</td>
  <td align="right">1.000000</td>
  <td align="right">+ 0.003000</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.995868</td>
  <td align="right">0.003004</td>
  <td align="right">0.998860</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">1</td>
  <td align="right">0.001039</td>
  <td align="right">0.062922</td>
  <td align="right">0.999447</td>
  <td align="right">0.062859</td>
  <td align="right">0.993123</td>
  <td align="right">0.002997</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.991754</td>
  <td align="right">5.275715</td>
  <td align="right">0.262476</td>
</tr>

<tr>
  <td align="center">1&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.981896</td>
  <td align="right">0.002999</td>
  <td align="right">0.984840</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">2</td>
  <td align="right">0.008211</td>
  <td align="right">0.125843</td>
  <td align="right">0.997789</td>
  <td align="right">0.125718</td>
  <td align="right">0.972886</td>
  <td align="right">0.002989</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.967552</td>
  <td align="right">2.605474</td>
  <td align="right">0.518508</td>
</tr>

<tr>
  <td align="center">2&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.955465</td>
  <td align="right">0.002988</td>
  <td align="right">0.958319</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">3</td>
  <td align="right">0.027155</td>
  <td align="right">0.188765</td>
  <td align="right">0.995025</td>
  <td align="right">0.188577</td>
  <td align="right">0.940420</td>
  <td align="right">0.002975</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.928937</td>
  <td align="right">1.701968</td>
  <td align="right">0.762083</td>
</tr>

<tr>
  <td align="center">3&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.917695</td>
  <td align="right">0.002971</td>
  <td align="right">0.920421</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">4</td>
  <td align="right">0.062586</td>
  <td align="right">0.251686</td>
  <td align="right">0.991155</td>
  <td align="right">0.251437</td>
  <td align="right">0.897462</td>
  <td align="right">0.002956</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.878252</td>
  <td align="right">1.241164</td>
  <td align="right">0.988001</td>
</tr>

<tr>
  <td align="center" colspan="13"><font size="+1"><b>&vellip;</b></td>
</tr>

<tr>
  <td align="center">95</td>
  <td align="right">6.769823</td>
  <td align="right">5.977546</td>
  <td align="right">- 3.988996</td>
  <td align="right">6.001390</td>
  <td align="right">0.0002917</td>
  <td align="right">- 0.021945</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.000057</td>
  <td align="right">0.000422</td>
  <td align="right">0.189466</td>
</tr>

<tr>
  <td align="center">95&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.0002844</td>
  <td align="right">- 0.021852</td>
  <td align="right">0.0002782</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">96</td>
  <td align="right">6.777942</td>
  <td align="right">6.040467</td>
  <td align="right">- 4.094580</td>
  <td align="right">6.065200</td>
  <td align="right">0.0002773</td>
  <td align="right">- 0.022473</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.000054</td>
  <td align="right">0.000405</td>
  <td align="right">0.185762</td>
</tr>

<tr>
  <td align="center">96&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.002705</td>
  <td align="right">- 0.022366</td>
  <td align="right">0.0002644</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">97</td>
  <td align="right">6.785828</td>
  <td align="right">6.103389</td>
  <td align="right">- 4.201270</td>
  <td align="right">6.129031</td>
  <td align="right">0.0002638</td>
  <td align="right">- 0.023006</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.000051</td>
  <td align="right">0.000389</td>
  <td align="right">0.182163</td>
</tr>

<tr>
  <td align="center">97&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.0002574</td>
  <td align="right">- 0.022884</td>
  <td align="right"> 0.0002515 </td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">98</td>
  <td align="right">6.793488</td>
  <td align="right">6.166310</td>
  <td align="right">- 4.309066</td>
  <td align="right">6.192881</td>
  <td align="right">0.0002511</td>
  <td align="right">- 0.023545</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.000048</td>
  <td align="right">0.000374</td>
  <td align="right">0.178666</td>
</tr>

<tr>
  <td align="center">98&frac12;</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.0002450</td>
  <td align="right">- 0.023408</td>
  <td align="right"> 0.0002393 </td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">99</td>
  <td align="right">6.800930</td>
  <td align="right">6.229232</td>
  <td align="right">- 4.417967</td>
  <td align="right">6.256752</td>
  <td align="right">0.0002391</td>
  <td align="right">- 0.024090</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="center">---</td>
  <td align="right">0.000045</td>
  <td align="right">0.000359</td>
  <td align="right">0.175267</td>
</tr>
</table>


----

[[File:DensityPert51Core02.png|325px|right|Density perturbation]]Hence, from the [[SSC/Perturbations#Continuity_Equation|linearized continuity equation]],

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

<tr>
  <td align="right">
<math>d</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- r_0 \frac{dx}{dr_0} - 3 x 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl\{
\biggl(\frac{4\pi}{45}\biggr)r_0^2  
- 3 \biggl[1 - \biggl(\frac{2\pi}{45}\biggr)r_0^2  \biggr]
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl\{
2\biggl(\frac{2\pi}{45}\biggr)r_0^2  
-3  + 3 \biggl(\frac{2\pi}{45}\biggr)r_0^2  
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl\{
\biggl(\frac{2\pi}{9}\biggr)r_0^2  
-3    
\biggr\} \, .
</math>
  </td>
</tr>
</table>
The solid black curve in the figure shown here, on the right, displays this analytically specified density perturbation, <math>d(m)</math>, for the case <math>\alpha_\mathrm{scale} = - 0.001</math>.

<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<div align="center"><b>Finite-Difference Representation of <math>d</math></b></div>
Our finite-difference representation of the mass-density at each radial shell in the equilibrium configuration (subscript "0") is,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[{\bar\rho}_{j-1/2} \biggr]_0</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[m_j - m_{j-1}\biggr]_0\biggl[\frac{4\pi}{3} (r_j^3 - r_{j-1}^3)  \biggr]_0^{-1} \, .
</math>
  </td>
</tr>
</table>
After perturbing the radial location of each shell &#8212; that is, after setting <math>r_j = (r_0)_j + [x \cdot r_0]_j</math> &#8212; the resulting finite-difference representation of the perturbed mass density of each shell is given by the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>{\bar\rho}_{j-1/2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[m_j - m_{j-1}\biggr]_0\biggl[\frac{4\pi}{3} (r_j^3 - r_{j-1}^3)  \biggr]^{-1} \, .
</math>
  </td>
</tr>
</table>
(Note that we retain a subscript "0" on the mass, <math>m_j</math>, because it serves as our Lagrangian identifier for each shell.)  The fractional density perturbation at each discrete shell is, then,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>
d_{j-1/2} \equiv \biggl[ \frac{\delta \rho}{\rho_0} \biggr]_{j-1/2}
=
\biggl\{ \frac{ {\bar\rho}_{j-1 / 2} }{ [{\bar\rho}_{j-1 / 2}]_0 } - 1 \biggr\}
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\frac{(r_j^3 - r_{j-1}^3)_0}{(r_j^3 - r_{j-1}^3)}  -1 
\, .
</math>
  </td>
</tr>
</table>
The red dots in the above "density perturbation" figure display how <math>d_j</math> varies with mass shell when <math>x</math> is specified by the parabolic dependence on radius.  The  dots lie virtually on top of the solid black curve in this figure, indicating that our finite-difference representation of the perturbed mass density matches well the analytically specified perturbed mass density.
</td></tr></table>

And, from [[SSC/Perturbations#Entropy_Conservation|linearized entropy conservation]],

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

<tr>
  <td align="right">
<math>p</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl\{
\biggl(\frac{2\pi}{9}\biggr)r_0^2  
-3    
\biggr\} \gamma_c 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl\{\biggl(\frac{2\pi}{9}\biggr)r_0^2  -3    \biggr\} \frac{6}{5}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \frac{dp}{dr_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl(\frac{8\pi}{15}\biggr)r_0  
\, .
</math>
  </td>
</tr>
</table>

Now, [[#Model_Amodel2|from below]] we find that,
<table border="0" cellpadding="8" align="center">

<tr>
  <td align="right"><math>\frac{P_0}{\rho_0} = \frac{P^*}{\rho^*}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl(1 + \frac{\xi^2}{3}\biggr)^{-3} \biggl(1 + \frac{\xi^2}{3}\biggr)^{5 / 2}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl(1 + \frac{\xi^2}{3}\biggr)^{-1 / 2} \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right"><math>g_0 = \frac{m}{(r^*)^2}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl( \frac{2\cdot 3}{\pi } \biggr)^{1/2} \biggl[ \xi^3 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] 
\biggl[ \biggl(\frac{3}{2\pi}\biggr)^{1 / 2}\xi \biggr]^{-2}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \biggl[ \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] 
\, .
</math>
  </td>
</tr>
</table>

Hence, the [[SSC/Perturbations#Euler_+_Poisson_Equations|linearized Euler + Poisson Equations]] expression gives,

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

<tr>
  <td align="right">
<math>\omega^2 r_0 x</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{P_0}{\rho_0} \frac{dp}{dr_0} 
- (4x + p)g_0 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl(1 + \frac{\xi^2}{3}\biggr)^{-1 / 2}  
\biggl(\frac{8\pi}{15}\biggr)r_0 
- \biggl\{
4\alpha_\mathrm{scale}\biggl[1 - \biggl(\frac{2\pi}{45}\biggr)r_0^2  \biggr] 
+ 
\alpha_\mathrm{scale} \biggl[\biggl(\frac{2\pi}{9}\biggr)r_0^2  -3    \biggr] \frac{6}{5}
\biggr\}
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \biggl[ \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \biggl(\frac{\omega^2}{\alpha_\mathrm{scale}}\biggr)\biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{3/2}  r_0 x</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(1 + \frac{\xi^2}{3}\biggr)  
\biggl(\frac{8\pi}{15}\biggr)r_0 
- \biggl\{
4 - \biggl(\frac{8\pi}{45}\biggr)r_0^2  
+ \biggl(\frac{12\pi}{45}\biggr)r_0^2  - \frac{18}{5} 
\biggr\}
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \xi  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(1 + \frac{\xi^2}{3}\biggr)  
\biggl(\frac{8\pi}{15}\biggr)r_0 
- \biggl[
\frac{2}{5}   
+ \biggl(\frac{4\pi}{45}\biggr)r_0^2  
\biggr]
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \xi  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(1 + \frac{\xi^2}{3}\biggr)  
\biggl(\frac{8\pi}{15}\biggr)\biggl(\frac{3}{2\pi}\biggr)^{1 / 2} \xi 
- \biggl[
\frac{2}{5}   
+ \biggl(\frac{4\pi}{45}\biggr)\biggl(\frac{3}{2\pi}\biggr) \xi^2  
\biggr]
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \xi  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(\frac{2^5\pi}{3 \cdot 5^2}\biggr)^{1 / 2} \xi 
- \frac{2}{5} \biggl( \frac{8\pi}{3 } \biggr)^{1/2} \xi  
+
\biggl(\frac{2^5\pi}{3^3\cdot 5^2}\biggr)^{1 / 2} \xi^3 
- \biggl( \frac{2^5\pi}{3^3\cdot 5^2 } \biggr)^{1/2} \xi^3  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>0</math>&nbsp; &nbsp; &nbsp; &nbsp; <font color="red">Exactly!</font>
  </td>
</tr>
</table>

====Step 6====

<font color="red"><b>STEP 6:</b></font> &nbsp; Given that when a proper solution has been obtained,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{d^2r}{dt^2}</math>
  </td>
  <td align="center">
<math>\rightarrow</math>
  </td>
  <td align="left">
<math>
- \omega^2 \cdot {\delta r}_i
\, ,
</math>
  </td>
</tr>
</table>
at each radial shell we can determine what the value of <math>\omega^2</math> would be as a result of our <math>{\delta r}_i</math> ''guess'' by rewriting the 

<div align="center">
<span id="PGE:Euler"><font color="#770000">'''Euler + Poisson Equations'''</font></span><br />
<math>\frac{d^2 r}{dt^2} = - 4\pi r^2 \frac{dP}{dM_r} - \frac{GM_r}{r^2} </math><br />
</div>
to read,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>+\biggl[ \frac{\omega_i^2}{G\rho_c} \biggr]</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\frac{1}{{\delta r}_i \cdot r_i^2} \biggl[ 4\pi (r_i^*)^4 \cdot \frac{dP^*}{dm^*} + m^*\biggr] \, .
</math>
  </td>
</tr>
</table>

<table border="1" align="center" cellpadding="5%" width="80%">
<tr><td align="left">

<div align="center">
<span id="PGE:Euler"><font color="#770000">'''Euler + Poisson Equations'''</font></span>
</div>

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

<tr>
  <td align="right">
<math>\biggl[ \frac{K_c^{1 / 2}}{G^{1 / 2} \rho_c^{2 / 5}} \biggr] \frac{d^2 r^*}{dt^2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 4\pi (r^*)^2 \biggl[ \frac{K_c^{1 / 2}}{G^{1 / 2} \rho_c^{2 / 5}} \biggr]^2 
\frac{dP^*}{dm} \biggl[ K_c \rho_c^{6 / 5} \biggr] \biggl[ \frac{G^{3 / 2}\rho_c^{1 / 5}}{K_c^{3 / 2}} \biggr]
- \frac{Gm}{(r^*)^2} \biggl[ \frac{K_c^{3 / 2}}{G^{3 / 2}\rho_c^{1 / 5}} \biggr]
\biggl[ \frac{G^{1 / 2} \rho_c^{2 / 5}}{K_c^{1 / 2}} \biggr]^2</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow~~~ \biggl[ \frac{K_c^{1 / 2}}{G^{1 / 2} \rho_c^{2 / 5}} \biggr] \frac{d^2 r^*}{dt^2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- 4\pi (r^*)^2 \biggl[ K_c^{1 / 2}G^{1 / 2}\rho_c^{3 / 5}\biggr] \frac{dP^*}{dm} 
- \frac{m}{(r^*)^2} \biggl[ K_c^{1 / 2}G^{1 / 2}\rho_c^{3 / 5} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow~~~ -\biggl[ \frac{1}{G \rho_c} \biggr] \frac{d^2 r^*}{dt^2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
4\pi (r^*)^2 \frac{dP^*}{dm} 
+ \frac{m^*}{(r^*)^2} 
</math>
  </td>
</tr>
</table>

----

<table align="center" border="0" cellpadding="5">
<tr>
  <td align="right">
<math>
\frac{d^2r}{dt^2}  
</math>
  </td>
  <td align="center">
<math>
\rightarrow
</math>
  </td>
  <td align="left">
<math>
\frac{d}{dt}\biggl[(i\omega) r_0 x~e^{i\omega t}\biggr] = - \omega^2 r_0 x~e^{i\omega t}
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
r^2 \frac{dP}{dm}  
</math>
  </td>
  <td align="center">
<math>
\rightarrow
</math>
  </td>
  <td align="left">
<math>
r_0^2 \biggl[1 + x~ e^{i\omega t}  \biggr]^2 \biggl\{\frac{dP_0}{dm} \biggl[1 + p~ e^{i\omega t}  \biggr]  + P_0~e^{i\omega t} \frac{dp}{dm}   \biggr\} \approx r_0^2 \frac{dP_0}{dm} \biggl[1 + (2x+p)~ e^{i\omega t}  \biggr] + P_0 r_0^2~e^{i\omega t} \frac{dp}{dm}
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>
\frac{Gm}{r^2}  
</math>
  </td>
  <td align="center">
<math>
\rightarrow
</math>
  </td>
  <td align="left">
<math>
\frac{Gm}{ r_0^2} \biggl[1 + x~ e^{i\omega t}  \biggr]^{-2} \approx \frac{Gm}{ r_0^2} 
\biggl[1 -2 x~ e^{i\omega t}  \biggr] \, .
</math>
  </td>
</tr>
</table>

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

After introducing a perturbation, we find that &hellip;

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

<tr>
  <td align="right">
RHS
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
4\pi \biggl[ r^* \biggr]^2 \frac{dP^*}{dm} 
+ m^*\biggl[ r^* \biggr]^{-2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
4\pi r_0^2\biggl[ 1 + \frac{\delta r}{r_0} \biggr]^2 \frac{d(P_0 + \delta P)}{dm} 
+ \frac{m^*}{r_0^2}\biggl[ 1 + \frac{\delta r}{r_0} \biggr]^{-2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
4\pi r_0^2\biggl[ 1 + 2x\biggr] \frac{d(P_0 + \delta P)}{dm} 
+ g_0\biggl[ 1 - 2x \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
4\pi r_0^2\biggl[ \frac{d(P_0 + \delta P)}{dm}\biggr]
+ g_0
+
4\pi r_0^2\biggl[ 2x\biggr] \frac{d(P_0 + \cancelto{\mathrm{small}}{\delta P})}{dm} 
-2x g_0 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
4\pi r_0^2\biggl[ \frac{d(P_0 + \delta P)}{dm}\biggr]
+ g_0
+ 2x
\biggl[ 4\pi r_0^2 \frac{dP_0}{dm} 
-g_0 \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>\biggl\{
4\pi r_0^2\biggl[ \frac{d(P_0 + \delta P)}{dm}\biggr]
+ g_0 \biggr\}
- 4xg_0
\, .
</math>
  </td>
</tr>
</table>

This "RHS" expression must be paired with &hellip;

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

<tr>
  <td align="right">
LHS
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[\frac{\omega_i^2}{G\rho_c}\biggr] \delta r_i
=
\biggl[\frac{\omega_i^2}{G\rho_c}\biggr] x r_0 \, .
</math>
  </td>
</tr>
</table>

The term inside the curly braces on the RHS is easy to evaluate inside our finite-difference scheme.  For our example parabolic displacement function, we expect this term to give,

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

<tr>
  <td align="right">
<math>\biggl\{
4\pi r_0^2\biggl[ \frac{d(P_0 + \delta P)}{dm}\biggr]
+ g_0 \biggr\}
</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\frac{P_0}{\rho_0} \cdot \frac{dp}{dr_0} - pg_0
=
4xg_0
\, .
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl(1 + \frac{\xi^2}{3}\biggr)^{-1 / 2}  
\biggl(\frac{8\pi}{15}\biggr)r_0 
- \biggl\{
\alpha_\mathrm{scale} \biggl[\biggl(\frac{2\pi}{9}\biggr)r_0^2  -3    \biggr] \frac{6}{5}
\biggr\}
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \biggl[ \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \alpha_\mathrm{scale} \biggl[ \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr]r_0^2  
\biggl( \frac{2^7\pi^3 }{3^3 \cdot 5^2 } \biggr)^{1/2}  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale}  
\biggl( \frac{2^5 \cdot 3^3\pi}{5^2 } \biggr)^{1/2}  \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2}  
+
\alpha_\mathrm{scale}\biggl(\frac{2^3\pi}{3\cdot 5}\biggr) \biggl(1 + \frac{\xi^2}{3}\biggr)^{-1 / 2}  r_0  
- 
\alpha_\mathrm{scale} r_0^2 
\biggl( \frac{2^7\pi^3}{3^3\cdot 5^2 } \biggr)^{1/2} \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2}  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggl\{
\biggl( \frac{2^5 \cdot 3^3\pi}{5^2 } \biggr)^{1/2}  \xi   
+
\biggl(\frac{2^5\pi}{3\cdot 5^2}\biggr)^{1 / 2} \biggl(1 + \frac{\xi^2}{3}\biggr) \xi   
- 
\biggl( \frac{2^5\pi}{3 \cdot 5^2 } \biggr)^{1/2} \xi^3 
\biggr\}  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggl\{
\biggl[ \biggl( \frac{2^5 \cdot 3^3\pi}{5^2 } \biggr)^{1/2}   
+
\biggl(\frac{2^5\pi}{3\cdot 5^2}\biggr)^{1 / 2}\biggr] \xi   
+
\biggl[ \biggl(\frac{2^5\pi}{3^3\cdot 5^2}\biggr)^{1 / 2}    
- 
\biggl( \frac{2^5\cdot 3^2\pi}{3^3 \cdot 5^2 } \biggr)^{1/2}\biggr] \xi^3 
\biggr\}  
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale} \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggl\{
\biggl(\frac{2^7\pi}{3}\biggr)^{1 / 2} \xi   
- 2\biggl(\frac{2^5\pi}{3^3\cdot 5^2}\biggr)^{1 / 2}  \xi^3 
\biggr\}  \, .
</math>
  </td>
</tr>
</table>
For comparison, note that this expression is identical to the expression for <math>4xg_0</math>, namely,

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

<tr>
  <td align="right"><math>4xg_0 </math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
4 \alpha_\mathrm{scale} \biggl[1 - \frac{\xi^2}{15}\biggr] 
\biggl( \frac{8\pi}{3 } \biggr)^{1/2} \biggl[ \xi \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\alpha_\mathrm{scale}\biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2}  
\biggl\{
\biggl( \frac{2^7\pi}{3 } \biggr)^{1/2}\xi - \biggl( \frac{2^7\pi}{3^3\cdot 5^2 } \biggr)^{1/2}\xi^3\biggr\} 
\, .
</math>
  </td>
</tr>
</table>

====Step 7====

<font color="red"><b>STEP 7:</b></font> &nbsp; Alternatively, from our [[SSC/Perturbations#Summary_Set_of_Linearized_Equations|summary set of linearized equations]], we expect &hellip;
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>
\omega^2 r_0 x
</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\frac{P_0}{\rho_0} \frac{dp}{dr_0} - (4x + p)g_0 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{P_0}{\rho_0} \cdot \frac{d(\delta P/P_0)}{dr_0} - \biggl(\frac{4\delta r}{r_0} + \frac{\delta P}{P_0}\biggr)\frac{Gm}{r_0^2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{P_0}{\rho_0} \biggl[ \frac{1}{P_0}\cdot \frac{d(\delta P)}{dr_0} - \frac{\delta P}{P_0^2} \frac{dP_0}{dr_0}\biggr]
+ 
\biggl[\biggl(1 - \frac{4\delta r}{r_0}\biggr) - \biggl(1 + \frac{\delta P}{P_0}\biggr)\biggr]\frac{Gm}{r_0^2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{\rho_0} \biggl[ \frac{d(\delta P)}{dr_0} - \frac{\delta P}{P_0} \frac{dP_0}{dr_0}\biggr]
+ 
\biggl(1 - \frac{4\delta r}{r_0}\biggr)\frac{Gm}{r_0^2} + \biggl(1 + \frac{\delta P}{P_0}\biggr) \frac{1}{\rho_0} \cdot \frac{dP_0}{dr_0}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{\rho_0} \biggl[ \frac{d(P_0 + \delta P)}{dr_0} \biggr]
+ 
\biggl[1 + \frac{\delta r}{r_0}\biggr]^{-4}\frac{Gm}{r_0^2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{\rho_0} \biggl[ \frac{dP}{dr_0} \biggr]
+ 
\frac{r_0^4}{r^4}\frac{Gm}{r_0^2} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{r_0^2}{r^4} \biggl[
4\pi r^4 \cdot \frac{dP}{dm}
+ 
Gm 
\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\Rightarrow ~~~ \omega^2
</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\frac{1}{\delta r \cdot r^2} \biggl[
4\pi r^4 \cdot \frac{dP}{dm}
+ 
Gm 
\biggr]\frac{r_0^2}{r^2} 
</math>
  </td>
</tr>
</table>

====Step 8====

<font color="red"><b>STEP 8:</b></font> &nbsp; Finally, we should ''guess'' a new eigenvector (then guess again, and again, and again &hellip;) until <math>\omega^2</math> settles down to have the same value at all radial locations.