> 
restart:with (plots): B:=1:M:=1:Gr:=0.5:Pr:=3:w:=0.02:blt:=5:Bi:=10:

> 
Eq1:=diff(f(eta),eta,eta,eta)(diff(f(eta),eta))^(2)+f(eta)*diff(f(eta),eta,eta)+B*H(eta)*(F(eta)diff(f(eta),eta))M*diff(f(eta),eta)+Gr*theta(eta)=0;


(1) 
> 
Eq2:=(1+Nr)*diff(theta(eta),eta,eta)+Pr*f(eta)*diff(theta(eta),eta)+(2/3)*H(eta)*B*(theta1(eta)theta(eta))=0;


(2) 
> 
Eq3:=H(eta)*F(eta)+H(eta)*diff(G(eta),eta)+G(eta)*diff(H(eta),eta)=0;


(3) 
> 
Eq4:=F(eta)^2+G(eta)*diff(F(eta),eta)+B*(F(eta)diff(f(eta),eta))=0;


(4) 
> 
Eq5:=G(eta)*diff(G(eta),eta)+B*(f(eta)+G(eta))=0;


(5) 
> 
Eq6:=G(eta)*diff(theta1(eta),eta)+l*B*(theta1(eta)theta(eta))=0;


(6) 
> 
bcs:=f(0)=0,(D(f))(0)=1,(D(theta))(0)=Bi*(1theta(0)),(D(f))(blt)=0,F(blt)=0,G(blt)=f(blt),H(eta)=w,theta(blt)=0,theta1(blt)=0;


(7) 

(8) 
> 
for k from 1 to 4 do p:=dsolve(eval({Eq1,Eq2,Eq3,Eq4,Eq5,Eq6,bcs},Nr=L[k]),[f(eta),F(eta),G(eta),H(eta),theta(eta),theta1(eta)],numeric,output=listprocedure);end do:

> 

> 

