In[]:=
ClearAll[Ω,h,z,deq2,deq3,Ω0,σ0,σa,λ,σ,r0]​​Needs["ErrorBarPlots`"]​​λ=2.1​​r0=0.3089​​Ω0=
r0
+1​​σ0=-2​​σa=0.04​​σ[z]=σ0+σa*z/(1+z)​​​​deq2=-2*(1+z)*h[z]*D[h[z],z]+3*h[z]^2-4*h[z]*Ω[z]+Ω[z]^2+2*(1+z)*h[z]*D[Ω[z],z]​​deq3=
-1
3
*(4*h[z]*Ω[z]-Ω[z]^2-2*(1+z)*h[z]*D[Ω[z],z])+σ0+σa*
z
1+z
*(2*h[z]*Ω[z]-Ω[z]^2)+λ​​​​solutiontrial=ParametricNDSolvedeq20,
-1
3
*(4*h[z]*Ω[z]-Ω[z]^2-2*(1+z)*h[z]*D[Ω[z],z])+x+y*
z
1+z
*(2*h[z]*Ω[z]-Ω[z]^2)+b==0,h[0]1,Ω[0]Ω0,{h,Ω},{z,0,4},{x,y,b}​​solmanipulate=h/.solutiontrial;​​​​Quiet[Manipulate[Plot[solmanipulate[x,y,b][z],{z,0,3}],{x,-5,5},{y,-5,5},{b,-5,5}]]​​​​hlambda=
0.3075*(1+z)^3+0.6911
​​ilambda=Plot[hlambda,{z,0,3},FrameTrue,AxesFalse,FrameLabel{"z","h ( z )"},PlotStyleRed,PlotRangeAll]​​​​
Out[]=
2.1
Out[]=
0.3089
Out[]=
1.55579
Out[]=
-2
Out[]=
0.04
Out[]=
-2+
0.04z
1+z
Out[]=
3
2
h[z]
-4h[z]Ω[z]+
2
Ω[z]
-2(1+z)h[z]
′
h
[z]+2(1+z)h[z]
′
Ω
[z]
Out[]=
2.1+-2+
0.04z
1+z
(2h[z]Ω[z]-
2
Ω[z]
)+
1
3
(-4h[z]Ω[z]+
2
Ω[z]
+2(1+z)h[z]
′
Ω
[z])
Out[]=
hParametricFunction
Expression: h
Parameters: {x,y,b}
,ΩParametricFunction
Expression: Ω
Parameters: {x,y,b}

Out[]=
​
x
y
b
0.5
1.0
1.5
2.0
2.5
3.0
-1.0
-0.5
0.5
1.0
Out[]=
0.6911+0.3075
3
(1+z)
Out[]=