In[]:=
ClearAll["Global`*"]
In[]:=
TOVeqs[rho0_]:=NDSolve[{​​kappa0*rho'[r]+(M[r]/(2r^2))*(1+kappa0*rho[r])*(1+4*Pi*r^3*kappa0*rho[r]^2/M[r])*(1-2M[r]/r)^(-1)0,​​M'[r]4Pi*rho[r]*r^2,​​rho[0.001]rho0,​​M[0.001]4Pi*rho0*(0.001)^3/3,​​kappa0==138.8},​​{rho,M},​​{r,0.001,20}];​​rhoTOV[r_,rho0_]:=rho[r]/.TOVeqs[rho0];​​MassTOV[r_,rho0_]:=M[r]/.TOVeqs[rho0];​​Rstar[rho0_]:=FindRoot[rhoTOV[r,rho0]0,{r,0.001,20}][[1,2]];​​Mstar[rho0_]:=MassTOV[Rstar[rho0],rho0][[1]];​​DataMstarRStar=Table[{1.47664*Rstar[10^i],Mstar[10^i]},{i,-4,-1,0.01}];​​StabTOV=ListLinePlot[DataMstarRStar,PlotStyle{Black},FrameTrue,FrameLabel{"
R
*
[km]","M/
M
⊙
"},LabelStyleMedium,PlotRange{{6,24},{-0.1,2.5}}]​​​​
In[]:=
ClearAll["Global`*"]
In[]:=
TOVeqs[rho0_]:=NDSolve[{​​kappa0*rho'[r]+(M[r]*rho[r]^(1/3)/(r^2))*(1+kappa0*rho[r]^(2/3))*(1+4*Pi*r^3*kappa0*rho[r]^(5/3)/M[r])*(1-2M[r]/r)^(-1)0,​​M'[r]4Pi*rho[r]*r^2,​​rho[0.001]rho0,​​M[0.001]4Pi*rho0*(0.001)^3/3,​​kappa0==49.6},​​{rho,M},​​{r,0.001,20}];​​TOVeqs[0.1]
In[]:=
rhoTOV[r_,rho0_]:=rho[r]/.TOVeqs[rho0];​​MassTOV[r_,rho0_]:=M[r]/.TOVeqs[rho0];​​Plot[{rhoTOV[x,0.01]},{x,0.001,20}]
Rstar[rho0_]:=FindRoot[rhoTOV[r,rho0]10^{-7},{r,0.001,14}][[1,2]];​​Mstar[rho0_]:=MassTOV[Rstar[rho0],rho0][[1]];​​
In[]:=
DataMstarRStar=Table[{1.47664*Rstar[10^i],Mstar[10^i]},{i,-4,-1,0.1}];​​StabTOV=ListLinePlot[DataMstarRStar,PlotStyle{Black,Thick},FrameTrue,FrameLabel{"
R
*
[km]","M/
M
⊙
"},LabelStyleMedium,PlotRange{{6,24},{-0.1,2.5}}]