PART ONE (TWO EXAMPLES): THE LAME-EMDEN ORDINARY DIFFERENTIAL EQUATION 
(ODE) FOR POLYTROPIC STELLAR INTERIORS
----------------------------------------------------------------------

     OVERVIEW:  The general Lane-Emden ODE (y is the mass density divided by 
     mass density at the center of the star, all raised to the power 1/n (n 
     is an aribrary positive integer); x is distance from the center of the
     star multiplied by a scaling factor):

                    y''(x) + 2y'(x)/x + y^n = 0

     Symbol "^" means above that y is raised to the power n.

     Boundary conditions:

     y(0) = 1   <== scaled mass density at center of star (x=0) is 1
     y'(0) = 0  <== rate of mass density change with radius is zero
                    at center of star

     The integer parameter n is defined from the relation

        GAMMA = 1+1/n, where Pressure = Constant * Density^GAMMA


     EXAMPLE #1 - N=1.5 POLYTROPE (I.E., GAMMA = 5/3):
     -------------------------------------------------

     In the runge.f ODE solver program, try the following:

     Domain of integration, etc.:
     NUM = 10000
     A = 0.00001  <== note: x yields a discontinuity at zero for this ODE
     B = 15
     SKIP = 10

     Initial boundary conditions:
     U(1) = 0   <== i.e., value of y'(x) at beginning
     U(2) = 1   <== value of y(x) at beginning

     Function declarations:
     F1 = -2*U1/X - U2*SQRT(ABS(U2))
     F2 = U1
     F3 = 0
     F4 = 0
        . 
        .
        .
     F10 = 0

     Meaningful values for the solution of y(x) are for x=0 extending
     to the first zero value for y(x) (maximum x value is about 3.65).

     (NOTE: This GAMMA=5/3 polytrope model also applies to a perfect 
     adiabatic atmosphere model for the Earth (but with different boundary 
     conditions than above)


     EXAMPLE #2 - N=3 POLYTROPE (I.E., GAMMA = 4/3):
     -----------------------------------------------

     Domain of integration, etc.:
     NUM = 10000
     A = 0.00001
     B = 15
     SKIP = 10

     Initial boundary conditions:
     U(1) = 0
     U(2) = 1

     Function declarations:
     F1 = -2*U1/X - U2*U2*U2
     F2 = U1
     F3 = 0
     F4 = 0
        .
        .
        .
     F10 = 0

     Meaningful values for the solution of y(x) are for x=0 extending
     to the first zero value for y(x) (maximum x value is about 6.90).




PART TWO (A THIRD EXAMPLE):  USING THE CODE TO SOLVE COUPLED FIRST-ORDER 
ODE'S:
------------------------------------------------------------------------

     Consider three containers of liquid or gas of fixed total mass and 
     with constant mass transfer rates between the three containers:

             |------------|            |------------|
             |            |     r1     |            |
             |     U1     | ---------> |     U2     |
             |            |            |            |
             |------------|            |------------|
                     ^                     |
                  r3 |                     | r2
                     |    |------------|   |
                     ---- |            | <--
                          |     U3     |
                          |            |
                          |------------|

     The math ("X" is time (in minutes), "U" is mass (in kg), and r1,r2,r3
     are the rate coefficients (units kg per minute)):
        
                 dU1                        dU2
                 --- = -r1*U1 + r3*U3       --- = -r2*U2 + r1*U1
                 dX                         dX
        
                              dU3
                              --- = -r3*U3 + r2*U2
                              dX
                         
     This example has the following rate coefficients and initial (X = 0)
     values:
                 r1 = 1
                 r2 = 1
                 r3 = 2
            
                 U1(0) = U2(0) = U3(0) = 1
            
     Try setting the domain of integration, etc. in the ODE solver 
     program as follows:

                 NUM = 10000
                 A = 0
                 B = 15
                 SKIP = 10

     In the function calls (make sure that F4 = F5 = ... = F10 = 0):

                 F1 = -U1 + 2*U3
                 F2 = -U2 + U1
                 F3 = -2*U3 + U2

     Final values (as time X tends toward infinity):

                 U1(inf) = 1.2
                 U2(inf) = 1.2
                 U3(inf) = 0.6
            
     NOTE: These final values can also be determined by applying the simple
     steady-state condition for mass (dU1/dt = dU2/dt = dU3/dt = 0 as time
     X --> infinity):

                 U1 + U2 + U3 = 3  <== Total mass is 3 kg
            
                 -U1 + 2*U3 = 0    <== These follow from the above
                 -U2 + U1 = 0      <== coupled differential equations
                 -2*U3 + U2 = 0    <== under steady-state condition

     It does not matter how the initial 3 kg of mass is distributed - the
     final mass values (X --> infinity) are 1.2, 1.2, and 0.6 kg, respectively.