load'runge_kutta.ijs'

NB. Series expansion of the Lane-Emden equation: 
NB. Usage: n LEps xi --> theta(xi), d theta/d xi
NB. http://www.scientificarts.com/laneemden/Links/laneemden_lnk_1.html
LEps=: 4 : 0
c6=. (x%15120)* (5 _8) p. x
c8=. (x%3265920)* (70 _183 122) p. x
c10=. (x%1796256000)* (3150 _10805 12642 _5032) p. x
c12=. (x%840647808000)* (138600 _574850 915935 _663166 183616) p. x
coeff=. 1,0,(%_6),0,(x%120),0,c6,0,c8,0,c10,0,c12 
th=. coeff p. y                NB. expansion for theta
thp=. (}.(i. 13)*coeff) p. y   NB. derivative of expansion
th,thp
)

NB. This defines the Lane-Emden diff eqn, x = n
Poly=: 4 : 0
n=. x
'x y w'=. y
dwdx=. -((2*w%x)+y^n)
w,dwdx
)

NB. To integrate the polytropic equation with n=1.5 from 0.1 to 0.3:
NB. (1.5&Poly) odeint 0.1 0.998335 _0.0332834; 0.3; 1e_9 0.001 1e_10

NB. Using the series to start, n=2:
NB. (2&Poly) odeint (0.02,(2 LEps 0.02)); 1.5; 1e_12 0.001 1e_10
NB. The zero (edge) of the n=3 polytrope is near 6.8968486:
NB. (3&Poly) odeint (0.02,(3 LEps 0.02)); 6.89684862; 1e_12 0.001 1e_10
NB. (4&Poly) odeint (0.04,(4 LEps 0.04)); 14.97154635; 1e_12 0.001 1e_10
NB. (4.8&Poly) odeint (0.04,(4 LEps 0.04)); 83.8128; 1e_12 0.001 1e_10

NB. Integration of Lane-Emden equation:  n LE_Run step
NB. E.g., the n=3 case:  3 LE_Run 0.05  &  to see the results, try
NB. ('pensize 3') plot (Re 0{"1 z);(Re 1{"1 z) ; plot (Re 0{"1 z);(Re 3{"1 z)
LE_Run=: 4 : 0
n=. x
xi=. step=. y
DEq=: n&Poly
z=: 0,(theta=. 1),0,1     NB. rows of z are (xi,theta,theta',rho)
x=. n LEps xi
rho=. (0{x)^n
z=: >z; (xi, x, rho)
while. theta > 0 do.
xi=. xi + step
y =. DEq odeint (}: {:z);xi;1e_12 0.001 1e_10
rho=. (1{y)^n
z=: z,(y,rho)
theta=. Re 1{y
end.
ypenu=. }: _2{z
xi=. xi - -:step
y =. DEq odeint ypenu;xi;1e_12 0.001 1e_10
rho=. (1{y)^n
last=. y,rho
)