Which $\gamma$ are you using? Also, you may consider the radius at $\sim10^4g/cm^3$, where you should have $^{56}$Fe nuclei (they have the highest binding energy per nucleon) plus electrons. The outer crust consists of a lattice of iron-56. At $\sim10^4g/cm^3$ the atoms are fully ionized owing to the high density. In any case, try to see the difference in radius between the radius at $10^4g/cm^3$ and, say, $10^6g/cm^3$: it should be almost the same! Related: TOV equation and the equation of state, solving numerically

@Quillo, I'm using all $\gamma$, from $\gamma \approx 1$ up to $\gamma = 2$ (even if non-realistic, to show the effects). Since I'm not using any value of $\kappa$, the output does not depend on the matter type or the micro-physics (except from the $\gamma$). I'll find the radius value in km only if I switch back from unitless variables ($\tilde{p}$ and $\tilde{r}$) to unitfull variables ($p$ and $r$), according to the description below equation (4) (please, read again!).

Cham, the point of my comment is that you should be able to find the radius at a given finite (small) density that, if converted to dimensional units, is close to the one of the Iron layer (i.e., the "cutoff" you are mentioning at the end of the question). Then I remember that for certain polytropic indexes, the solution is unbounded, maybe it's worth to try also with gamma smaller than 1 or bigger than 2: do you know what happens in those cases? (for example, the ideal Fermi gas should have $\gamma=5/3<1$... if I remember correctly, it is now 20 years since I last solved the TOV numerically).

@Quillo, $\gamma = 5/3$ is certainly larger than 1! I'm not really interested in $\gamma < 1$ and $\gamma > 2$ which are physically very weird (my code shows all sorts of weird singularities with these values). I don't think the cutoff is related in any way to the iron layer. On the contrary, the cutoff is defined at very low pressure and density, which isn't in the core of the star.

Can you plot log p on the y-axis. It's impossible to see what's going on when the density should vary over many orders of magnitude.

Hauahaha sorry Cham, silly mistake! The rest is however a reasonable piece of advice. Also, I agree with ProfRob that a linear-log plot would be great. The cutoff is arbitrary, you can justify it only on the basis of physical arguments and the iron layer is a possibility that is often used. When you have the log plot it would be amazing to see it for two different numerical methods (as you probably know you can select the method in the DSolve function).

The log(p) plot is strange, it should go to minus infinity for the pressure approaching zero at the surface. Why do you use Method -> "StiffnessSwitching"? Did you try other methods? I think it's a numerical problem. My runge kutta 4 integrator for the TOV worked well. I was cutting at the iron layer because the last part of the star is completely irrelevant and tabulated realistic equatiosn of state never go so low in density anyway. However, the polytrope was working well with RK4.

Also,please make sure you are showing us an example with $\gamma >1.2$ because the radius is infinite for smaller values

@ProfRob, I updated my question with two plots for $\gamma = \frac{4}{3}$ (which is larger than your $\gamma = \frac{6}{5} = 1.20$ ($n = 5$). From this plot, I don't see how I could find a finite radius.

About a theoretical proof of infinite radius for $𝑛>5$, I found the following in my notes: For a newtonian polytropic fluid, we can get the potential energy:\begin{equation} 𝑈=−\, \frac{3}{5−𝑛} \, \frac{𝐺𝑀^2}{𝑅}.\end{equation} So this expression implies $𝑛<5$ to get $𝑈<0$ (gravity is attractive). It gives $𝑈>0$ for all $𝑛>5$, unless $R\rightarrow\infty$. So indeed $\gamma=1+\tfrac{1}{𝑛} = \tfrac{6}{5}=1.20$ is "special". I'm not sure this is also valid for the relativistic (non-newtonian) regime but it makes sense.

@ProfRob, how do we know that the infinite radius for $𝑛 > 5$ is also valid for the relativistic regime? The previous argument is from Newtonian's gravity only.

@Cham I don't know. Besides, that isn't your problem - because you've now checked that, right?

You have still not extended the x-axis of the logarithmic plot far enough for us to see what is happening as the density becomes small, even for the Newtonian case. Why is the "normalised" central pressure 0.5?

@ProfRob, the pressure at 0.5 is just one example. I can change that value from 0 to 10 (or even more), in the units chosen (for which $G \equiv 1$ and $\kappa \equiv 1$. The shape of the pressure curve is "typical".