The following code illustrates the principles of BVP solution.

It is simplified to be clear enough and thus it does not exploit the usually very high precision of the pseudospectral method

solveTOV[pCentral_Real,OptionsPattern[{γ->2.7,dim->9 ,init->{rMax->0.1,m[0]->0.01,Φ[0]->-0.01,p[0]->0.01}}]]:=Block[
{γEOS=OptionValue[γ],n=OptionValue[dim],initGuess=OptionValue[init], 
rT,sl,residualΦ,residualP,residualM,bdryConditions,chebRoots,rPts,eqns,unknowns,unknowns0,solution},

rT[n_Integer,r_]:=ChebyshevT[n,-1+2r/rMax]; 
chebRoots[n_Integer]:=Table[  Cos[  (2k+1) Pi/2/n],{k,n-1,0,-1}]//N;
rPts ={r-> rMax(#+1)/2}&/@chebRoots[n+1];

sl = {p->Sum[ p[k] rT[k,r],{k,0,n}],
Φ->Sum[Φ[k] rT[k,r],{k,0,n}],
m->Sum[m[k] rT[k,r],{k,0,n}]r^3
};
  eos = ρ-> HeavisideTheta[p]( Abs[p]^(1/γEOS)+p/(γEOS-1));

residualΦ =Hold[ D[Φ,r]-(m+4 Pi r^3 p)/r/(r-2m)  ]  /.eos/.sl//ReleaseHold;
residualP =Hold[ D[p,r] +(ρ+p)D[Φ,r]  ] /.eos/.sl//ReleaseHold;
residualM =Hold[ D[m,r] - 4Pi r^2 ρ ]  /.eos/.sl//ReleaseHold;

bdryConditions = Join[
{p}/.sl/.r->rMax,
{Exp[2Φ]-(1-2m/r)}/.sl/.r->rMax,
{p-pCentral}/.sl/.r->0
];

eqns = {residualM/.rPts[[2;;]],residualΦ/.rPts,residualP/.rPts[[;;-2]]}~Join~bdryConditions//Flatten;

unknowns = {rMax}~Join~Table[{Φ[k],p[k],m[k]},{k,0,n}]//Flatten;
unknowns0 = Transpose[ {unknowns,unknowns/.initGuess/.{rMax->0.1,m[0]->0.01,Φ[0]->-0.01,p[0]->0.01}/.{m[k_]->0,p[k_]->0,Φ[k_]->0}}];

solution = FindRoot[ eqns, unknowns0,MaxIterations->200];

Join[ solution , {M->(m/.sl/.r->rMax/.solution),p->(p/.sl/.solution),m->(m/.sl/.solution),Φ->(Φ/.sl/.solution)}]
]

With[ {pC=1.0},
tmp = solveTOV[pC,dim->16,init->{rMax->0.5,Φ[0]->-0.5}];
Print[Take[tmp,7]];
GraphicsColumn[{
Plot[m/.tmp,{r,0,rMax/.tmp},PlotLabel->m,PlotTheme->"Detailed",PlotLegends->None],
Plot[p/.tmp,{r,0,rMax/.tmp},PlotLabel->p,PlotTheme->"Detailed",PlotLegends->None],
Plot[Φ/.tmp,{r,0,rMax/.tmp},PlotLabel->Φ,PlotTheme->"Detailed",PlotLegends->None]}]
]

prevSol={rMax->0.1,m[0]-> 0.001,Φ[0]-> -0.001,p[0]-> 0.001};
Table[ {rMax,M}/. (solveTOV[pC,dim->8,γ->2.7,init->prevSol]),{pC,10^Range[-6,+0.5,0.125]}];
ListLinePlot[%]