load'/your_path/runge_kutta.ijs'

NB. 12-term 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=: dyad define 
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
th, (}.(i. 13)*coeff) p. y   NB. & derivative of expansion
)

NB. This defines the Lane-Emden diff eqn. 
NB. (x = the polytropic index n)
Poly=: dyad define
'xi theta w'=. y
dwdx=. -(2*w%xi)+ theta^x
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 to surface. 
NB. Usage:  n LE_Run step-size
NB. The result is stored in the global variable "z"
LE_Run=: dyad define
n=. x
xi=. step=. y
DEq=: n&Poly
z=: 0,(theta=. 1),0,1     NB. rows of z are (xi,theta,theta',rho)
yy=. n LEps xi
rho=. n^~ {.yy
z=: >z; xi, yy, rho 
params=. 1e_12, 0.001, 1e_10  NB. error, trial step, minimum step
while. theta > 0 do.
  xi=. xi+ step
  yy=. DEq odeint (}: {:z);xi; params
  rho=. (1{yy)^n
  z=: z,(yy,rho)
  theta=. Re 1{yy
end.
yult=: }: _1{z    NB. last step (}: trims off density)
ypenu=: }: _2{z   NB. penultimate step
NB. linear interpolation for xi such that theta(xi)=0
step0=. Re (1{ypenu)* -%/ }:yult-ypenu  
NB. past 0 we can get small imaginary result, hence "Re"
xi=. (xi-step)+ 0.997*step0  NB. 0.997: stop short of 0
yy=. DEq odeint ypenu;xi; params
rho=. (1{yy)^n
z=: Re (}:z),(last=. yy,rho) NB. gard against complex type   
last
)

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)  for theta(xi), or
NB. 'pensize' plot (Re 0{"1 z);(Re 3{"1 z)  for density