Neutron star modelling with polytrope

I want to model the density profile of a neutron star with 1.4 $$M_\odot$$ with a polytrope $$p=\kappa\rho^{\gamma}$$ with $$\gamma=1+\frac{1}{n}$$. I started with $$n=\frac{1}{2}$$ and solved the Lane-Emden equation for the Newtonian case with Mathematica. Then I put the resulting numerical solution for $$\rho(r)$$ into the structure equations for stars (with TOV equation for $$\frac{\text{d}p}{\text{d}r}$$) and solved them numerically.

For the boundary conditions I used $$R=13\,$$km and the conditions from the Lane-Emden equation: $$\kappa=\rho_c^{1-\frac{1}{n}}\frac{4\pi G}{n+1}\frac{1}{a^2}$$, $$a=\frac{\xi_0}{R}$$, $$\rho_c=\left(-\frac{3}{\xi_0}\partial_z\theta(\xi_0)\right)^{-1}\rho_{\text{avg}}$$, $$\rho_{\text{avg}}=\frac{3M}{4\pi R^3}$$. Where $$\xi_0$$ is the Root of $$\theta(z)$$.

Now when I put $$M=1.4\,M_{\odot}$$ the resulting mass after solving the TOV equation is $$\approx 0.9\,M_{\odot}$$. Only if I go up to around $$3.9\,M_{\odot}$$ I get the desired mass. Also changing the polytrope index $$n$$ doesn't really help.

Somehow I think, I have the wrong approach by mixing the Newtonian and full GR (TOV equation) case. Is there a way to get a realistic density profile with a polytrope and $$M=1.4\,M_{\odot}$$ by using the TOV equation?

The equation to look at is ($$2.8$$), which reduces to the regular (Newtonian) Lane-Emden form when the parameter $$\sigma \equiv \dfrac{P}{\rho c^2}$$ goes to zero.