γ = 2.5

Here right-hand side of differential equations is defined

ρ[p_Real] := Abs[p]^(1/γ)
μ[p_Real] := ρ[p]+p/(γ-1)
dpdr[r_Real,p_Real,m_Real] := If[r>0,-( μ[p]+p) (m+4 Pi r^3 p)/(r(r-2m)),0]

eqs = {
p'[r] ==dpdr[r,p[r],m[r]],
m'[r] ==4 Pi r^2  μ[p[r]]
}

We solve it as initial-value ODE problem and ignore ϕ(r)

sol = NDSolve[ eqs && {p[0]==1,m[0]==0}&& WhenEvent[p[r]<0 ,"StopIntegration"] &&
WhenEvent[ r-2m[r]<0,"StopIntegration"] ,{p[r],m[r]},{r,0,10}]

The results are plotted

Plot[1-2m[r]/r/.sol, {r,0,sol[[1,1,2,0]]["Domain"][[1,2]] },PlotTheme->"Detailed"]
Plot[p[r]/.sol, {r,0,sol[[1,1,2,0]]["Domain"][[1,2]] },PlotTheme->"Detailed"]

To plot mass-radius relation we need following function

rMaxM[p0_Real]:= Block[{sol,rMax},
sol = NDSolve[ eqs && {p[0]==p0,m[0]==0}&& WhenEvent[p[r]<0 ,"StopIntegration"] &&
WhenEvent[ r-2m[r]<0,"StopIntegration"] ,{p[r],m[r]},{r,0,10}];
rMax= sol[[1,1,2,0]]["Domain"][[1,2]];
{rMax,sol[[1,2,2,0]][rMax]}
]
(*rMaxM[1.]*)

ParametricPlot[ rMaxM[Exp[logpC]],{logpC,-15,15},PlotPoints->30,MaxRecursion->1,PlotRange->Full,AspectRatio->2/3,PlotTheme->"Detailed",FrameLabel->{r,m}]