@Preben Alsholm 

 Thank you for your help. I will try using your suggestion.

@Preben Alsholm 

restart: with(plots):blt:=15.5:A:=1:M:=1:lambda:=1:fw:=1:rhos:=8933:rhof:=997.1:rhob:=14918.11:rhob2:=20939.1:kn:=0.613:kf:=401:cps:=385:cpf:=4179:Pr:=6.2:
Eq1 := (1/(1-s)^(2.5))*diff(f(eta),eta,eta,eta)-(1-s+s*(rhos/rhof))*(diff(f(eta),eta)^(2)-f(eta)*diff(f(eta),eta,eta)+A*(diff(f(eta),eta)+(1/2)*eta*diff(f(eta),eta,eta)))-M*diff(f(eta),eta)+(1-s+s*((rhob/rhob2)))*lambda*theta(eta)=0;
d / d / d \\
----- |----- |----- f(eta)|| /
deta \ deta \ deta // |
---------------------------- - (1 + 7.958981045 s) |
2.5 \
(1 - s)

2
/ d \ / d / d \\
|----- f(eta)| - f(eta) |----- |----- f(eta)||
\ deta / \ deta \ deta //

\
/ d \ 1 / d / d \\|
+ |----- f(eta)| + - eta |----- |----- f(eta)|||
\ deta / 2 \ deta \ deta ///

/ d \
- |----- f(eta)| + (1 - 0.2875476978 s) theta(eta) = 0
\ deta /
Eq2 :=(1/Pr)*((kn/kf)/(1-s+s*(rhos*cps/rhof*cpf)))*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-diff(f(eta),eta)*theta(eta)-A*(2*theta(eta)+0.5*eta*diff(theta(eta),eta))=0;
/ d / d \\
0.0002465610168 |----- |----- theta(eta)||
\ deta \ deta //
------------------------------------------
7
1 + 1.441423799 10 s

/ d \ / d \
+ f(eta) |----- theta(eta)| - |----- f(eta)| theta(eta)
\ deta / \ deta /

/ d \
- 2 theta(eta) - 0.5 eta |----- theta(eta)| = 0
\ deta /
bcs1 := f(0) = fw, (D(f))(0) = 1, (D(f))(blt) = 0, theta(0) = 1, theta(blt) = 0;
f(0) = 1, D(f)(0) = 1, D(f)(15.5) = 0, theta(0) = 1,

theta(15.5) = 0
L := [0,0.05,0.1,0.2];
[0, 0.05, 0.1, 0.2]
for k to 4 do R := dsolve(eval({Eq1, Eq2, bcs1}, s = L[k]), [f(eta),theta(eta)], numeric, output = listprocedure);Y || k := rhs(R[3]); end do:
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

@Preben Alsholm