In[]:=
p[ϵ[r_]]:=K*
2
ϵ[r]
​​K=
8
10
;
In[]:=
eq1=λ'[r]==
1
4*r(1+2*α*R[r]+α*r*R'[r])
*2*(1+2*α*R[r])*(1-Exp[2*λ[r]])+
r²*Exp[2*λ[r]]
3
(16*Pi*(2*ϵ[r]+3*p[ϵ[r]])+2*R[r]+3*α*R[r]^2)+
2*α*r*R'[r]
1+2*α*R[r]
*2*(1+2*α*R[r])*(1-Exp[2*λ[r]])+
r²*Exp[2*λ[r]]
3
(16*Pi*ϵ[r]+R[r]+3*α*R[r]^2)+4*α*r*R'[r];
In[]:=
eq2=ψ'[r]==
1
4*r(1+2*α*R[r]+α*r*R'[r])
*(r²*Exp[2*λ[r]](16*Pi*p[ϵ[r]]-α*R[r]^2)+2*(1+2*α*R[r])*(Exp[2*λ[r]]-1)-8*α*r*R'[r]);
In[]:=
eq3=R''[r]==
Exp[2*λ[r]]
6*α
*(8*Pi*(3*p[ϵ[r]]-ϵ[r])+R[r])+R'[r]*λ'[r]-ψ'[r]-
2
r
;
In[]:=
eq4=D[p[ϵ[r]],r]==-(ϵ[r]+p[ϵ[r]])ψ'[r];
In[]:=
r0=
-20
10
;​​rg=1.477*1000(*m*);​​Rc=-(4.7)*
-3
10
*
-2
rg
;​​α=1*
-2
rg
;
In[]:=
sol1=NDSolve[{eq1,eq2,eq3,eq4,λ[r0]==0,ψ[r0]==0,R[r0]==Rc,R'[0]==0,ϵ[r0]=(2.0*7.4263)*
-10
10
},{λ,ψ,R,ϵ},{r,r0,20}]
NDSolve
:Equation or list of equations expected instead of 1.48526×
-9
10
in the first argument 
′
λ
[r]
2(1-Power[2])(1+9.16788×
-7
10
R[1])+
1
3
1r²(1)+
22r1(1)
1+Times[2]
4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])
,
′
ψ
[r]
1
1
,5,
′
R
[0]0,1.48526×
-9
10
.
NDSolve
:Equation or list of equations expected instead of 1.48526×
-9
10
in the first argument 
′
λ
[r]
2(1-Power[2])(1+9.16788×
-7
10
R[1])+
1
3
1r²(1)+
22r1(1)
1+Times[2]
4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])
,
′
ψ
[r]
1
1
,5,
′
R
[0]0,1.48526×
-9
10
.
Out[]=
NDSolve
′
λ
[r]2(1-
2λ[r]

)(1+9.16788×
-7
10
R[r])+
1
3
2λ[r]

r²2R[r]+1.37518×
-6
10
2
R[r]
+16π2ϵ[r]+300000000
2
ϵ[r]
+9.16788×
-7
10
r
′
R
[r]2(1-
2λ[r]

)(1+9.16788×
-7
10
R[r])+
1
3
2λ[r]

r²R[r]+1.37518×
-6
10
2
R[r]
+16πϵ[r]+1.83358×
-6
10
r
′
R
[r](1+9.16788×
-7
10
R[r])(4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])),
′
ψ
[r]2(-1+
2λ[r]

)(1+9.16788×
-7
10
R[r])+
2λ[r]

r²-4.58394×
-7
10
2
R[r]
+1600000000π
2
ϵ[r]
-3.66715×
-6
10
r
′
R
[r](4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])),
′′
R
[r]363588.
2λ[r]

R[r]+8π-ϵ[r]+300000000
2
ϵ[r]
+
′
R
[r]-
2
r
+
′
λ
[r]-
′
ψ
[r],200000000ϵ[r]
′
ϵ
[r]-ϵ[r]-100000000
2
ϵ[r]

′
ψ
[r],λ
1
100000000000000000000
0,ψ
1
100000000000000000000
0,R
1
100000000000000000000
-2.15445×
-9
10
,
′
R
[0]0,1.48526×
-9
10
,{λ,ψ,R,ϵ},r,
1
100000000000000000000
,20
In[]:=
sol1
Out[]=
NDSolve
′
λ
[r]2(1-
2λ[r]

)(1+9.16788×
-7
10
R[r])+
1
3
2λ[r]

r²2R[r]+1.37518×
-6
10
2
R[r]
+16π2ϵ[r]+300000000
2
ϵ[r]
+9.16788×
-7
10
r
′
R
[r]2(1-
2λ[r]

)(1+9.16788×
-7
10
R[r])+
1
3
2λ[r]

r²R[r]+1.37518×
-6
10
2
R[r]
+16πϵ[r]+1.83358×
-6
10
r
′
R
[r](1+9.16788×
-7
10
R[r])(4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])),
′
ψ
[r]2(-1+
2λ[r]

)(1+9.16788×
-7
10
R[r])+
2λ[r]

r²-4.58394×
-7
10
2
R[r]
+1600000000π
2
ϵ[r]
-3.66715×
-6
10
r
′
R
[r](4r(1+9.16788×
-7
10
R[r]+4.58394×
-7
10
r
′
R
[r])),
′′
R
[r]363588.
2λ[r]

R[r]+8π-ϵ[r]+300000000
2
ϵ[r]
+
′
R
[r]-
2
r
+
′
λ
[r]-
′
ψ
[r],200000000ϵ[r]
′
ϵ
[r]-ϵ[r]-100000000
2
ϵ[r]

′
ψ
[r],λ
1
100000000000000000000
0,ψ
1
100000000000000000000
0,R
1
100000000000000000000
-2.15445×
-9
10
,
′
R
[0]0,1.48526×
-9
10
,{λ,ψ,R,ϵ},r,
1
100000000000000000000
,20