In[]:=
f[r_]:=1-
2M
2
r
3/2
(
2
r
+
2
q
)
+
8πP
2
r
3
​​
In[]:=
Solve[f[r]==0,M]
Out[]=
M
3/2
(
2
q
+
2
r
)
(3+8Pπ
2
r
)
6
2
r

In[]:=
mass=
3/2
(
2
q
+
2
r
)
(3+8Pπ
2
r
)
6
2
r
;
In[]:=
rh=r/.Solve[f[r]==0,r,Reals];
In[]:=
temp=
D[f[r],r]
4π
/.Mmass//Simplify
Out[]=
-2
2
q
+
2
r
+8Pπ
4
r
4πr(
2
q
+
2
r
)
In[]:=
q=1;​​M=2;​​Plot[temp/.rrh,{P,0,0.3},PlotStyleRed,AxesOrigin{0,0}]
Out[]=
0.05
0.10
0.15
0.20
0.25
0.30
-1.0
-0.5
0.5
1.0