## 195 Reputation

12 years, 244 days

Amir

## What is the problem?...

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?

## Thanks...

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

## physical value...

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?

Dear Preben

## @Carl Love    Thanks Carl and ...

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));

## A Problem...

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));

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.

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)

Amir

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

## Finding the error...

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.

## A solution...

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

Amir

## A change in code...

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

Thanks Preben

## a Problem with algortihm...

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

## Singularity...

Thanks Peben and aslo Thanks to Carl.

Thanks

## a litte problem...

This method doesnt work when I have

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

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