In[]:=
eq1:=D[f[eta],{eta,3}]+(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][10]0,theta[10]0,phi[10]0};Pr=7;G[r]=5;G[m]=3;N[b]=0.1;N[t]=0.1;Le=1;gamma=0.2;tau=0.1;solForAlpha[alpha_]:=NDSolve[{eq1,eq2,eq3,ics,bcs}/.alphaalpha,{f,theta,phi},{eta,0,5},MaxStepSize->0.1,MaxSteps->10^3,Method{"StiffnessSwitching",Method{"ExplicitRungeKutta",Automatic}}];alphaValues={0°,30°,60°,90°};plotsf=Table[Plot[Evaluate[{f'[eta]}/.solForAlpha[alpha]],{eta,0,5},PlotStyleSwitch[alpha,0°,Brown,30°,Red,60°,Green,90°,Black],PlotLegends{StringJoin["α = ",ToString[alpha]]},AxesLabel{η,f'},PlotLabel"Velocity profile f'(η)"],{alpha,alphaValues}];Show[Flatten[plotsf],PlotRangeAll]
Out[]=
