Dear collagues

Hi,

I write a code to solve a system of ODE. It solve the ODES in a wide range of parameters but as I decrease NBT below 0.5, it doesnt converge. I do my best but I couldn't find the answer. Would you please help me? Thank you

Here is my code and it should be run for 0.1<NBT<10. the value of NBT is input directly in res1.

restart:

EPSILONE:=1000:

Digits:=15:

a[mu1]:=5.45:

b[mu1]:=108.2:

a[k1]:=1.292:

b[k1]:=-11.99:

rhop:=4175:

rhobf:=998.2:

mu1[bf]:=9.93/10000:

k1[bf]:=0.597:

rhost(eta):=1-phi(eta)+phi(eta)*rhop/rhobf;

k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);

eq1:=(diff(u(eta), eta))*a[mu1]*(diff(phi(eta), eta))+2*(diff(u(eta), eta))*b[mu1]*phi(eta)*(diff(phi(eta), eta))+((diff(u(eta), eta))+(diff(u(eta), eta))*a[mu1]*phi(eta)+(diff(u(eta), eta))*b[mu1]*phi(eta)^2)/(eta+EPSILONE)+diff(u(eta), eta, eta)+(diff(u(eta), eta, eta))*a[mu1]*phi(eta)+(diff(u(eta), eta, eta))*b[mu1]*phi(eta)^2+1-phi(eta)+phi(eta)*rhop/rhobf:

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

eq3:=diff(phi(eta),eta)+phi(eta)/(N_bt*(1+gama1*T(eta))^2)*diff(T(eta),eta):

eq2:=subs(phi(0)=phi0,eq2):

eq3:=subs(phi(0)=phi0,eq3):

eval(dsolve({eq3,phi(1)=phiv},phi(eta)),(T)(1)=1):

Phi:=normal(combine(%)):

Teq:=isolate(eval(eq2,Phi),diff(T(eta),eta)):

ueq1:=eval(eq1,Phi)=0:

ueq2:=subs(Teq,ueq1):

lambda:=0;

Ha:=0;

N_bt:=cc*NBT+(1-cc)*0.8;

kratio:=k1[p]/k1[bf]:

GUESS:=[T(eta) =0.0001*eta, u(eta) =0.1*eta, phi(eta) = 0.3*(eta-1)^4];

res1 := dsolve(subs(NBT=0.48,gama1=0.2,phiv=0.06,{eq1,eq2,eq3,u(0)=lambda*D(u)(0),D(u)(1)=0,T(0)=0,phi(1)=phiv,T(1)=1}), numeric,method=bvp[midrich],maxmesh=4000,approxsoln=GUESS, output=listprocedure,continuation=cc):

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

masst:=evalf(int((1-G0(eta)+G0(eta)*rhop/rhobf)*G1(eta), eta = 0..1));

heatt:=(1+a[k1]*G0(0)+b[k1]*G0(0)^2)*G2(0);

plots:-odeplot(res1,[[eta,T(eta)]],0..1,legend=[T],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);

plots:-odeplot(res1,[[eta,u(eta)]],0..1,legend=[u],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);

plots:-odeplot(res1,[[eta,phi(eta)]],0..1,legend=[phi],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);

>

>

Thank you

**Amir**