In[]:=
f[r_,P_]:=1-+
2M
2
r
3/2
(+)
2
r
2
q
8πP
2
r
3
In[]:=
q=1;M=3.5;
In[]:=
Clear[fPlot];fPlot[P_]:=D[f[r,P],r]/(4*Pi)/.Last@NSolve[f[r,P]==0&&r>0,r,Reals];
In[]:=
Plot[fPlot[P],{P,0,0.3},PlotRange->{0,0.5},PlotStyle->Green,AxesOrigin->{0,0}]
Out[]=