In[]:=
ClearAll["`*"];Remove["`*"]
In[]:=
m=2;phi=Rationalize[0.1,0];rhoF:=Rationalize[997.1,0];CpF:=4179;kF:=Rationalize[0.613,0];rhoS:=5200;CpS:=670;kS:=6;kNF:=kF*(kS+(m-1)*kF-(m-1)*phi*(kF-kS))/(kS+(m-1)*kF+phi*(kF-kS));alpha=30Degree;M=1;Pr=Rationalize[6.2,0];
In[]:=
sys={x1'[eta]==x2[eta],x2'[eta]==x3[eta],x3'[eta]==-((1-phi)^(5/2)((1-phi)rhoF+phi*rhoS)(1/2*x3[eta](x1[eta]Sin[alpha]+eta*Cos[alpha])-M(x2[eta]-1))),x4'[eta]==x5[eta],x5'[eta]==-(kF*Pr((1-phi)rhoF*CpF+phi*rhoS*CpS)(1/2(x1[eta]Sin[alpha]+eta*Cos[alpha])x5[eta]/kNF))};
In[]:=
pfun=With[{eps=0},ParametricNDSolveValue[{sys,x1[eps]==eps,x2[0]==eps,x3[0]==a,x4[eps]==1,x5[eps]==b},{x1[eta],x2[eta],x3[eta],x4[eta],x5[eta]},{eta,eps,1},{a,b},WorkingPrecision->20,Method"StiffnessSwitching"]];
In[]:=
Table[Plot[Evaluate[pfun[1/2,1/10][[j]]],{eta,0,1},PlotLabelsSubscript[x,j],WorkingPrecision->20],{j,5}](*For:a=1/2andb=1/10*)
Out[]=
,
,
,
,
Plot[Evaluate[pfun[1/2,1/10][[1]]],{eta,0,1},PlotLabelsSubscript[x,1][eta]](*For:a=1/2andb=1/10*)
Out[]=
Plot[Evaluate[pfun[1/2,1/10][[2]]],{eta,0,1},PlotLabelsSubscript[x,2][eta]](*For:a=1/2andb=1/10*)
Out[]=
Plot[Evaluate[pfun[1/2,1/10][[3]]],{eta,0,1},PlotLabelsSubscript[x,3][eta]](*For:a=1/2andb=1/10*)
Out[]=
Plot[Evaluate[pfun[1/2,1/10][[4]]],{eta,0,1},PlotLabelsSubscript[x,4][eta],WorkingPrecision->20](*For:a=1/2andb=1/10*)
Out[]=
Plot[Evaluate[pfun[1/2,1/10][[5]]],{eta,0,1},PlotLabelsSubscript[x,5][eta],WorkingPrecision->20](*For:a=1/2andb=1/10*)
Out[]=