Hello, I think your equations are quite if not too complicated. And your code does not run on my system (version 7.0).

So I changed the 1st part a bit. Could you compare your results to these? I get a negative rho (density?), which should be impossible.

kappa0 = 138.8

tov[rho0_] := Module[{},
  NDSolve[{kappa0*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\ \(rho[
          r]\)\) + (M[r]/(2 r^2))*(1 + kappa0*rho[r])*(1 + 
          4*Pi*r^3*kappa0*rho[r]^2/M[r])*(1 - 2 M[r]/r)^(-1) == 0, \!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\ \(M[r]\)\) == 4 Pi*rho[r]*r^2,
      rho[0.001] == rho0, M[0.001] == 4 Pi*rho0*(0.001)^3/3}, {rho, 
     M}, {r, 0.001, 20}] // Flatten]

rhoTOV[r_, rho0_] := rho[r] /. tov[rho0];
MassTOV[r_, rho0_] := M[r] /. tov[rho0];
Plot[{rhoTOV[x, 0.01]}, {x, 0.001, 20}, PlotRange -> All]