load'/your_path/runge_kutta.ijs'

NB. Series expansion of the isothermal sphere equation: 
NB. Usage: Iso_ser z --> w(z), dw/dz
Iso_ser=: monad define
coeff=. 0,0,(%6),0,(%120),0,(%,1890) 
w=. coeff p. y                NB. expansion for w(z)
wp=. (}.(i. 7)*coeff) p. y   NB. derivative of expansion
w,wp
)

NB. Equation for isothermal sphere
Iso_Sph=: monad define
'z yy w'=. y
dwdz=. (^-yy)-(2*w%z)
w,dwdz
)

NB. To integrate the isothermal equation from 0.1 to 5:
NB. get starting w and dw/dz: 
NB. 'w wp'=. Iso_ser z=. 0.1
NB. Iso_Sph odeint (z,w,wp) ; 5; 1e_9 0.001 1e_10

NB. Directly incorporating the starting series, 
NB. Iso_Sph odeint (0.02,(Iso_ser 0.02)); 10; 1e_12 0.001 1e_10

NB. Integration of Isothermal equation: Iso_Run step
NB. Iso_Run 0.05 200  &  to see the results, try
NB. ('pensize 3') plot (0{"1 zz);(1{"1 zz) ; plot (0{"1 zz);(3{"1 zz)
DEq=: Iso_Sph
Iso_Run=: monad define
'step max'=. y
z=. step
zz=: 0,0,0,1,0,0     NB. rows of zz are (z,w,w',^-w)
'w wp'=. Iso_ser z
rho=. ^-w
p=. sV=. 0
zz=: >zz; (z,w,wp,rho,p,sV)
ind=. _1
while. (ind=. >:ind)<max do. 
z=. z + step
V=. 'z w wp' =. DEq odeint (_3}. {:zz);z;1e_12 0.001 1e_10
rho=. ^-w
p=. (z^4)*(wp^2)*rho
sV=.% (z*wp)^3
zz=: zz,(V,rho,p,sV)
end.
last=. V,rho,p,sV
)

NB. Iso_Run 0.05 500
NB. rsq=. (10^. 1<. 2%}. *: 0{"1 zz)  
NB. rsq is the log of 1/z^2 but capped at 1 -- compare to iostherm:
NB. 'pensize 3' plot ( }. 0{"1 zz);>(10^. }. 3{"1 zz);rsq
NB. Iso_Run 0.05 5000
NB. rsq=. (10^. 1<. 2%}. *: 0{"1 zz) 
NB. 'pensize 3' plot ( }. 0{"1 zz);>(10^. }. 3{"1 zz);rsq

NB.  ww=. 2}. w=. 1{"1 zz
NB. 'pensize 3' plot ( 150{. rr);>(5%~ 150{. ww);(150{. rho) 
NB. 'pensize 3' plot ( 120{. rr);>(120{. ww);(120{. rho)
NB. 'pensize 3' plot (10^. 2000{. rr);>(10^. 2000{. ww);(10^. 2000{. rho)

NB.   Iso_Run 0.05 8000
NB. 400.05 11.2375 0.00497983 1.31715e_5