@Preben Alsholm 

Thanks.

It may help but does not improve the solution very much. I do not know what is the problem?

Why maple cannot give the solution for lower values of NBT? Is there a singularity or ...?

In addition, I observed some thing that may help. I printed the values of a[1] and a[2] and observed that these values reduced to about 1e-12. but, the solution continues and when these values are in the range of 0.01, get the results. It means that the solution procedure does not get the best results. Is there any chance to add a condition that when the values of a[1] and a[2] are about 1e-12, the solution get the resutls?

@Preben Alsholm 

Thanks preben. It improve my algorithm. However, as you mentioned we cannot go further NBT=0.44!

I dont know why we cannot get the results for NBt=0.1 and lower. (for some other cases, for example phiavg=0.02 and Ha=5 I get the results upto NBT=0.2)

What is wrong with the equations?

Is there any chance to get the results for lower values of NBT? for example change of variables or other numerical methods?

It should be mentioned that most of this code was written with your previous helps and I am very grateful for your helps.

Sincerely yours

Amir

@Preben Alsholm 

Thanks Preben.

Actually, phi(1) should be (physically) positive.

Is there any chance to force the code that dont give the negative values of phi?

what about the convergence region? is there a chance the wide the region of convergence?

Thanks for your attention in advance

@Preben Alsholm 

Dear Preben

Many thanks for your help.

@Carl Love 

Thanks Carl and Preben.

Actually, I have two problem

1. first, how canI change back fy(y) to f(x)? I want to plot f(x)-x after solving the odes.

2. this code still has a problem. When you chose inf=15 or 30 it doesnt converge.

In the following, I put the final code. thanks for your attentions in advance

restart:

# parametrs

MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:

dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:

k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:

phi:=0.0:
inf:=15: Pr:=7: lambda:=0.1:


with(plots):

PDEtools:-dchange({f(x)=fy(y),T(x)=Ty(y),h(x)=hy(y),x=y*inf},[dqu1=0,dqu2=0,dqu3=0],[fy,Ty,hy,y]);
sysy:=solve(%,{diff(fy(y),y,y,y),diff(Ty(y),y,y),diff(hy(y),y)});
#The boundary conditions:
bcsy:={Ty(0)=1,Ty(1)=0,fy(0)=0,D(fy)(0)=inf*lambda,D(fy)(1)=0,hy(1)=0};
indets(sysy,name);



pppe:=dsolve( {sysy[1],sysy[2],sysy[3],bcsy[1],bcsy[2],bcsy[3],bcsy[4],bcsy[5],bcsy[6]}, numeric );  

print(odeplot(pppe,[diff(f(x),x),x],0..inf,color=black,numpoints=400));

@Preben Alsholm 

Thanks Carl and Preben.

Actually, I have two problem

1. first, how canI change back fy(y) to f(x)? I want to plot f(x)-x after solving the odes.

2. this code still has a problem. When you chose inf=15 or 30 it doesnt converge.

In the following, I put the final code. thanks for your attentions in advance

restart:

# parametrs

MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:

dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:

k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:

phi:=0.0:
inf:=15: Pr:=7: lambda:=0.1:


with(plots):

PDEtools:-dchange({f(x)=fy(y),T(x)=Ty(y),h(x)=hy(y),x=y*inf},[dqu1=0,dqu2=0,dqu3=0],[fy,Ty,hy,y]);
sysy:=solve(%,{diff(fy(y),y,y,y),diff(Ty(y),y,y),diff(hy(y),y)});
#The boundary conditions:
bcsy:={Ty(0)=1,Ty(1)=0,fy(0)=0,D(fy)(0)=inf*lambda,D(fy)(1)=0,hy(1)=0};
indets(sysy,name);



pppe:=dsolve( {sysy[1],sysy[2],sysy[3],bcsy[1],bcsy[2],bcsy[3],bcsy[4],bcsy[5],bcsy[6]}, numeric );  

print(odeplot(pppe,[diff(f(x),x),x],0..inf,color=black,numpoints=400));

@Preben Alsholm

Dear Preben

when the code is run for each values parameters (NBT and Lambda), I want to store the value of Nub computed for each parameter in B[jj] (see the following):

G0,G1,G2:=op(subs(RES,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..1)))/(Phiavg*rhop+(1-Phiavg)*rhobf):

phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) :
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1)):
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..1))):
Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2))):
B[jj]:=Nub;

and also I want to plot this

odeplot(RES,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);

But I couldnt add them correctly.

would you please help me to add the following commands to the code.

I them with as follows

with(plots):

PP:=proc(Nbt,lambda,s,fi0) option remember; local u1,res,a,INT0,INT10 ;
global G0,G1,G2;
userinfo(2,procname,printf("N[bt] = %a lambda = %a sigma = %a phi0 = %a\n",Nbt,lambda,s,fi0));
if not type([Nbt,lambda,s,fi0],list(numeric)) then return 'procname(_passed)' end if;
u1:=u1sol(Nbt,lambda,s,fi0);
RES2(parameters=[Nbt,lambda,s,fi0,u1]); #Setting parameters
G0,G1,G2:=op(subs(RES2,[u(eta),phi(eta),T(eta)])):
end proc:

for jj from 1 by 1 to 54 do
print(jj);



RES2:=dsolve({ueq2,Teq,Phi,T(0)=0,u(0)=lambda*u1,D(u)(0)=u1},numeric, output=listprocedure,parameters=[N[bt],lambda,sigma,phi0,u1],maxfun=0);
#u1:=u1sol(A[jj](1),A[jj](2),A[jj](3),A[jj](4));
#RES2(parameters=[A[jj](1),A[jj](2),A[jj](3),A[jj](4),u1]); #Setting parameters
 PP(A[jj](1),A[jj](2),A[jj](3),A[jj](4));
#RES1:=dsolve({subs({N[bt]=A[jj](1),lambda=A[jj](2),sigma=A[jj](3),phi0=A[jj](4)},ueq2),subs({N[bt]=A[jj](1),lambda=A[jj](2),sigma=A[jj](3),phi0=A[jj](4)},Teq),subs({N[bt]=A[jj](1),lambda=A[jj](2),sigma=A[jj](3),phi0=A[jj](4)},eq3=0),T(0)=0,u(0)=A[jj](2)*D(u)(0),u(1)=-A[jj](2)*D(u)(1), phi(0)=A[jj](4)},output=listprocedure,numeric):

ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..1)))/(Phiavg*rhop+(1-Phiavg)*rhobf):
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) :
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1)):
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..1))):
Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2))):
B[jj]:=Nub;

odeplot(RES2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
forget(PP);
end do;

however, It gives wrong results. (I coudnt find the reason)

Thanks for your attention in advance

Amir

@Preben Alsholm 

Many Thanks Preben for your kind help.  You improve my code very well.

@Preben Alsholm 

Thanks Preben,

Actually, I dont know why this problem is tricky?

what's wrong with the problem which is very sensitive to these parameters?

This code can be very useful to me to find the solution of many type of novel problems which is not investigated before. So, I really eager to find a stronger algorithm which can easily solve these problems.

It is worth mentioning that your answer to my previous questions was very helpful, Especially using of continious funtion help me very well to find the solution of many problems easily.

Thanks for your attention in advance

Thanks Preben.

I am so sorry for the output in the code. I didnt know that.

Actually, I think the code must be changed in a way that the "complex value encountered" vanished.

For the case above, I added the Digits:=15; and the code worked for NBT=0.5, However, it doesnt work for NBT=0.4 !!!

Why the code is very sensitive to these options?

For each run, I must change many parameters in order to run the code.

I would be most grateful if you could guide me how to change the variables to get the good results?

I do not have any idea why the code is very  sensitive.

the final version of the code is :

restart:
Digits:=15;
Phiavg:=0.06;
lambda:=0.05;
Ha:=0;
NBT:=0.5;
Nr:=500;
#N[bt]:=cc*NBT+(1-cc)*4; ## cc between 0 and 1
N[bt]:=cc*NBT+(1-cc^2)*2;



eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(sigma-Nr*(phi(eta)-phi1[w])-(1-phi(eta))*T(eta)-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta)-1/(k(eta)/k1[w]);
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):
gama1:=0.00:
a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=1:
phi1[w]:=phi0:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);
Q:=proc(pp2,fi0) option remember; local res,F0,F1,F2,a,INT0,INT10,B;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(sigma=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}), numeric,output=listprocedure,initmesh=10, continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
INT0:=evalf(Int((abs(F0(eta)),eta=0..1)));
INT10:=evalf(Int(abs(F0(eta))*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf),eta=0..1));
#a[1]:=evalf(Int((F0(eta),eta=0..1)));
a[2]:=(INT10/INT0-Phiavg)/Phiavg*100; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0) Q(_passed)[2] end proc;
#Q(116,0.0041);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[130,0.01]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[43.55,0.39]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5.65,0.000036]);
tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[12,0.003]); # khoob ba 1
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[20,0.01]);
sigma:=tempe[2](1);
phi0:=tempe[2](2);
with(plots):

res2 := dsolve({eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}, numeric,output=listprocedure,continuation=cc);
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) ;
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet))):

Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));
(rhs(res2(0.0000000000001)[3])-rhs(res2(0)[3]))/0.0000000000001;
sigma;
NBT;

I want to run this code for 0.2<NBT<10 and lambda=0,0.05,0.1,0.2

Thanks for your attention in advance

Amir

@Preben Alsholm 

Thanks preben

for a[1] and a[2], for the relative measure I change

a[1]:=(B-10000*pp2)/10000/pp2 - 1  ; #relative
a[2]:=(INT10/INT0-Phiavg)/Phiavg -1 ; #relative

to

a[1]:=(B-10000*pp2)/10000/pp2; #relative
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative

I Think it is ok now.

Thanks

@Preben Alsholm 

Thanks Preben

Actually I think there is a problem with my algorithm.

I think it is very weak.

for example. consider the following code (not singular)

restart:

Phiavg:=0.02;
lambda:=0.0;
Ha:=0;
N[bt]:=0.4:
                              0.02
                               0.
                               0



eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(1-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*((1)*rho(eta)*c(eta)*u(eta)/(p2*10000)+( (a[k1]+2*b[k1]*phi(eta))/(1+a[k1]*phi1[w]+b[k1]*phi1[w]^2)*diff(phi(eta),eta)*diff(T(eta),eta) ));
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
          /  d   /  d         \\   mu1[w]
          |----- |----- u(eta)|| + -------
          \ deta \ deta       //   mu(eta)

                      /  d           \ /  d         \
               mu_phi |----- phi(eta)| |----- u(eta)|
                      \ deta         / \ deta       /
             + --------------------------------------
                              mu(eta)                
                                /      /                      
                                |      |                      
/  d   /  d         \\     1    |      |rho(eta) c(eta) u(eta)
|----- |----- T(eta)|| + ------ |k1[w] |----------------------
\ deta \ deta       //   k(eta) |      |       10000 p2       
                                \      \                      

                                /  d           \ /  d         \\\
     (a[k1] + 2 b[k1] phi(eta)) |----- phi(eta)| |----- T(eta)|||
                                \ deta         / \ deta       /||
   + ----------------------------------------------------------||
                                                  2            ||
                 1 + a[k1] phi1[w] + b[k1] phi1[w]             //
                                             /  d         \
                        2.500000000 phi(eta) |----- T(eta)|
     /  d           \                        \ deta       /
     |----- phi(eta)| - -----------------------------------
     \ deta         /                             2        
                                (1 - gama1 T(eta))         




mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):
gama1:=0.00:
a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=1:
#phi(0):=1:
#u(0):=0:
phi1[w]:=phi0:

mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);

/  d   /  d         \\   /                                    
|----- |----- u(eta)|| + \0.0009930000000 + 0.03883623000 phi0
\ deta \ deta       //                                        

                      2\//               
   + 0.5301627000 phi0 / \0.0009930000000

                                                   2\   
   + 0.03883623000 phi(eta) + 0.5301627000 phi(eta) / +

  /                                       /  d           \ /  d  
  |(0.03883623000 + 1.060325400 phi(eta)) |----- phi(eta)| |-----
  \                                       \ deta         / \ deta

         \\//                                        
   u(eta)|| \0.0009930000000 + 0.03883623000 phi(eta)
         //                                          

                          2\
   + 0.5301627000 phi(eta) /
                                                     /      
                                                     |      
   /  d   /  d         \\              1             |      
   |----- |----- T(eta)|| + ------------------------ |(0.597
   \ deta \ deta       //   0.597 + 4.45959 phi(eta) \      

                      /
                      |
                      |
      + 4.45959 phi0) |
                      \

     /             6                        6\       
     \-1.1752324 10  phi(eta) + 4.1744724 10 / u(eta)
     ------------------------------------------------
                         10000 p2                    

             /  d           \ /  d         \\\
        7.47 |----- phi(eta)| |----- T(eta)|||
             \ deta         / \ deta       /||
      + ------------------------------------||
                   1 + 7.47 phi0            //
     /  d           \                        /  d         \
     |----- phi(eta)| - 2.500000000 phi(eta) |----- T(eta)|
     \ deta         /                        \ deta       /
#method= bvp[midrich]
#A somewhat speedier version uses the fact that you really need only compute 2 integrals not 3, since one of the integrals can be written as a linear combination of the other 2:
Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10,B;
global Q1,Q2;
print(pp2,fi0);
#lambda/(phi(1)*rhop/rhobf+(1-phi(1)))
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve(subs(p2=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,D(u)(1)=0,u(0)=lambda*D(u)(0),phi(0)=phi0,D(T)(0)=1,T(0)=0}), numeric,output=listprocedure):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
INT0:=evalf(Int((F0(eta),eta=0..1)));
INT10:=evalf(Int(F0(eta)*F1(eta),eta=0..1));
B:=(-cbf*rhobf+cp*rhop)*INT10+ rhobf*cbf*INT0;
a[1]:=B-10000*pp2;
a[2]:=INT10/INT0-Phiavg;
Q1(_passed):=a[1];
Q2(_passed):=a[2];
if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if
end proc;
#The result agrees very well with the fsolve result.
#Now I did use a better initial point. But if I start with the same as in fsolve I get the same result in just about 2 minutes, i.e. more than 20 times as fast as fsolve:

Q1:=proc(pp2,fi0) Q[1](_passed) end proc;
Q2:=proc(pp2,fi0) Q[2](_passed) end proc;
evalf[5](exp(+1/N[bt]));
Optimization:-LSSolve([Q1,Q2],initialpoint=[115,0.0041]);



proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;

tempe:=%:

p2:=convert(tempe[2](1),float,18);

phi0:=convert(tempe[2](2),float,18);


                      116.436227708060514
                     0.00619399727967237176


with(plots):

res2 := dsolve(subs(p2=convert(tempe[2](1),float,12),phi0=convert(tempe[2](2),float,12),{eq1=0,eq2=0,eq3=0,D(u)(1)=0,u(0)=lambda*D(u)(0),phi(0)=phi0,D(T)(0)=1,T(0)=0}), numeric,output=listprocedure):
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(G0(eta)*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(G0(eta)*G1(eta),eta=0..1))) / evalf((Int(G0(eta),eta=0..1))) ;
TTb:=evalf(Int(G0(eta)*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(G0(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);



res2(0)[4];
  

h:=(4/TTb);
                   HFloat(8.903283428616747)
((1+a[k1]*G1(0)+b[k1]*G1(0)^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2));
                   HFloat(0.9102741949531519)
(4/TTb)*(((1+a[k1]*G1(0)+b[k1]*G1(0)^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));

It is run and get the accurate results. however, When I change the value of N[bt]:=0.2 It doenst converge.

Can I strenght the algorithm to give me the results for different values of N[bt]?

Amir

@Preben Alsholm 

Thanks Peben and aslo Thanks to Carl.

your comments about method= bvp[midrich] is very helpful.

Thanks

This method doesnt work when I have

D(T)(0)=1 (replace with T(0)=0).

Please see the following code:

restart:
eq1:=diff(u(eta),eta,eta)-dp+Gr[T]*T(eta)-Gr[phi]*phi(eta);
eq2:=diff(eta*(1+phi(eta)+phi(eta)^2)*T(eta),eta);
eq3:=Db*diff(phi(eta),eta)+Dt*diff(T(eta),eta);

bcs1 := phi(0)=phi[w];
bcs2 := D(T)(0)=1;
bcs3 := u(h)=0,D(u)(0)=0;

sol:=dsolve({eq3=0,bcs1},phi(eta));
res1:=eval(eval(sol,{bcs1}),{bcs2});
eq:=simplify(eval(eq2=0,res1));
res_T:=dsolve({eq,bcs2},T(eta));
res_phi:=eval(res1,res_T);
odetest([res_T,res_phi],[sys_ode]);
eq_u:=simplify(eval(eval(eq1=0,res_T),res_phi));
dsolve({eq_u=0,bcs3},u(eta));
res_u:=eval(%,{bcs3});



the above algorithm cannot remove the T(0) from sol and res1. also, it gave an error that i couldnt fix it

Thanks for your attention in advance

Amir