f[r_]:=1-
2M
r
+
2
Q
2
r
+
8πP
2
r
3
-
a
4
Q
20
6
r
-
c
1+3ω
r
;​​Tinv2=
2/3
2
a
4/3
π
4
Q
9
1/3
3
7/3
V
-
2
Q
3V
+
1
2
1/3
6
2/3
π
1/3
V
+P
1/3
6
π
1/3
V
+
-
2
3
+2ω
2
c
1
3
-ω
3
π
-
2
3
-ω
V
ω;​​ri=r/.Solve[-5a
4
Q
-4
4
r
(8πP-3
2
Q
+2
2
r
)-6c
5-3ω
r
ω(5+3ω)==0,r];​​rh=r/.Solve[f[r]==0,r];​​TI2=Tinv2/.V
4π
3
r
3
/.rri;​​Plot[TI2/.{a0.8,c0.2,ω-0.5,Q0.5},{P,0,1},PlotStyleBlack,AxesOrigin{0,0}]​​Plot
D[f[r],r]
4π
/.rrh/.{a0.8,c0.2,ω-0.5,Q0.5,M10},{P,0,5},PlotRange->{0,10},PlotStyleGreen,AxesOrigin{0,0}​​