In[]:=
ClearAll["`*"]
(*Thisdefinesthecoefficientsfora4(5)explicitRunge–KuttapairofFehlbergthatwaspopularinthe1960s*)
In[]:=
Fehlbergamat={​​{1/4},​​{3/32,9/32},​​{1932/2197,-7200/2197,7296/2197},{439/216,-8,3680/513,-845/4104},{-8/27,2,-3544/2565,1859/4104,-11/40}};​​Fehlbergbvec={25/216,0,1408/2565,2197/4104,-1/5,0};​​Fehlbergcvec={1/4,3/8,12/13,1,1/2};​​Fehlbergevec={-1/360,0,128/4275,2197/75240,-1/50,-2/55};​​FehlbergCoefficients[4,p_]:=​​N[{Fehlbergamat,Fehlbergbvec,Fehlbergcvec,Fehlbergevec},p];
In[]:=
phi=1/2;rhoF=1;rhoS=3/2;alpha=30Degree;M=1/10;kNF=1;kF=1/2;Pr=7/10;CpF=1;CpS=12/10;
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))}//Simplify;
In[]:=
pfun=ParametricNDSolveValue[{sys,x1[0]==0,x2[0]==0,x3[0]==a,x4[0]==1,x5[0]==b},{x1[eta],x2[eta],x3[eta],x4[eta],x5[eta]},{eta,0,1},{a,b},Method->{"ExplicitRungeKutta","Coefficients"FehlbergCoefficients,"DifferenceOrder"4,"EmbeddedDifferenceOrder"5,"StiffnessTest"False}];
Table[Plot[Evaluate[(pfun[a,b]/.a->1/.b->1)[[j]]],{eta,0,1},PlotLabelsSubscript[x,j]],{j,5}](*For:a=1andb=1*)
Out[]=

,
,
,
,

X1=1;Plot[Evaluate[Table[pfun[a,1][[X1]],{a,-1,1,1/4}]],{eta,0,1},PlotRangeAll](*Solutionforx1[eta],a=(-1,1),b=1*)
Out[]=
In[]:=
X4=4;Plot[Evaluate[Table[pfun[1,b][[X4]],{b,-1,1,1/4}]],{eta,0,1},PlotRangeAll](*Solutionforx4[eta],a=1,b=(-1,1)*)
Out[]=
X5=5;Plot[Evaluate[Table[pfun[a,b][[X5]],{a,-1,1,1/4},{b,-1,1,1/4}]],{eta,0,1},PlotRangeAll](*Solutionforx5[eta],a=(-1,1),b=(-1,1)*)
Out[]=