eq1:=(beta+1)*D[f[eta],{eta,3}]+beta*((1/2)*eta*Cos[alpha]*D[f[eta],{eta,2}]+(1/2)*Sin[alpha]*f[eta]*D[f[eta],{eta,2}]+G[r]*theta[eta]+G[m]*phi[eta])0;eq2:=D[theta[eta],{eta,2}]+(1/2)*Pr*(Cos[alpha]*eta*D[theta[eta],eta])+(1/2)*Pr*(Sin[alpha]*f[eta]*D[theta[eta],eta])+N[b]*D[theta[eta],eta]*D[phi[eta],eta]+N[t]*D[theta[eta],eta]^20;eq3:=D[phi[eta],{eta,2}]+(1/2)*Le*(Cos[alpha]*eta*D[phi[eta],eta])+(1/2)*Le*(Sin[alpha]*f[eta]*D[phi[eta],eta])+N[t]*D[theta[eta],{eta,2}]/N[b]0;ics:={f[0]0,Derivative[1][f][0]gamma*Derivative[2][f][0],theta[0]1+tau*Derivative[1][theta][0],phi[0]1};bcs:={Derivative[1][f][7]0,theta[7]0,phi[7]0};Pr=7;G[r]=5;G[m]=3;N[b]=0.1;N[t]=0.1;Le=1;gamma=0.1;tau=0.01;beta=1;
solForAlphaNN[alpha_]:=NDSolve[{eq1,eq2,eq3,ics,bcs}/.alphaalpha,{f,theta,phi},{eta,0,7},MaxStepSize->0.5,MaxSteps->10^3,Method{"StiffnessSwitching",Method{"ExplicitRungeKutta",Automatic}}];alphaValues={0°,30°,60°,90°};
plotsfNN=Table[Plot[Evaluate[{f'[eta]}/.solForAlphaNN[alpha]],{eta,0,7},PlotStyleSwitch[alpha,0°,Blue,30°,Red,60°,Green,90°,Black],PlotLegends{StringJoin["α = ",ToString[alpha]]},AxesLabel{η,f'},PlotLabel"Velocity profile f'(η)"],{alpha,alphaValues}];Show[Flatten[plotsfNN],PlotRangeAll]plotsthetaNN=Table[Plot[Evaluate[{theta[eta]}/.solForAlphaNN[alpha]],{eta,0,2.5},PlotStyleSwitch[alpha,0°,Blue,30°,Red,60°,Green,90°,Black],PlotLegends{StringJoin["α = ",ToString[alpha]]},AxesLabel{η,θ},PlotLabel"Temperature profile θ(η)"],{alpha,alphaValues}];Show[Flatten[plotsthetaNN],PlotRangeAll]plotsphiNN=Table[Plot[Evaluate[{phi[eta]}/.solForAlphaNN[alpha]],{eta,0,5},PlotStyleSwitch[alpha,0°,Brown,30°,Red,60°,Green,90°,Black],PlotLegends{StringJoin["α = ",ToString[alpha]]},AxesLabel{η,ϕ},PlotLabel"Concentration profile ϕ(η)"],{alpha,alphaValues}];Show[Flatten[plotsphiNN],PlotRangeAll]
Out[]=
Out[]=
Out[]=
