Dear collegues

I wrote the following code

restart:

Digits := 15;

a[k]:=0;

b[k]:=7.47;

a[mu]:=39.11;

b[mu]:=533.9;

mu[bf]:=9.93/10000;

k[bf]:=0.597;

ro[p]:=3880 ;

ro[bf]:= 998.2;

c[p]:= 773;

c[bf]:= 4182;

#mu[bf]:=1;

Gr[phi]:=0; Gr[T]:=0;

#dp:=0.1;

Ree:=1;

Pr:=1;

Nbt:=cc*NBTT+(1-cc^2)*6;

#######################

slip:=0.1; ####

NBTT:=2; ####

lambda:=0.1; ####

phi_avg:=0.02; ####

#######################

eq1:=diff( (1+a[mu]*phi(eta)+b[mu]*phi(eta)^2)*diff(u(eta),eta),eta)+dp/mu[bf]+Gr[T]*T(eta)-Gr[phi]*phi(eta);

eq2:=diff((1+a[k]*phi(eta)+b[k]*phi(eta)^2)*diff(T(eta),eta),eta)+lambda*T(eta)/k[bf];

eq3:=diff(phi(eta),eta)+1/Nbt*diff(T(eta),eta);

Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10;

global Q1,Q2;

print(pp2,fi0);

if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:

res := dsolve({subs(dp=pp2,eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=fi0}, numeric,output=listprocedure,continuation=cc);

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

a[1]:=evalf(Int(F0(eta),eta=0..1))-Ree*Pr;;

a[2]:=INT10/INT0-phi_avg;

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;

Q1:=proc(pp2,fi0) Q[1](_passed) end proc;

Q2:=proc(pp2,fi0) Q[2](_passed) end proc;

Optimization:-LSSolve([Q1,Q2],initialpoint=[0.3,0.0007]);

se:=%[2];

res2 := dsolve({subs(dp=se[1],eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=se[2]}, numeric,output=listprocedure,continuation=cc);

G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):

TTb:=evalf(Int(G0(eta)*G2(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1))/evalf(Int(G0(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1));

with(plots):

odeplot(res2,[[eta,phi(eta)/phi_avg]],0..1);

odeplot(res2,[[eta,T(eta)/TTb]],0..1);

odeplot(res2,[[eta,u(eta)/(Ree*Pr)]],0..1);

res2(1);

Nuu:=(1/TTb);

1/((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2));

(1/TTb)*(((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2)));

>

I want to run the code for the value of **NBTT in the range of 0.2 to 10. **this code gave the results in the range of 4-10 easily. So, I used the continuation which improve the range of the results between 2-10. However, I coudnt gave the results when 0.2<**NBTT**<2. Would you please help me in this situation.

Also, It is to be said that the values of phi should be positive. in some ranges, I can see that phi(1) is negative. Can I place a condition in which the values phi restricted to be positive.

Thanks for your attentions in advance

**Amir**