Dear Collagues

I wrote a code to solve a system of ode numerically. When the equations have been solved, I want to find the value of MTE which is in integral form. I check and see that the value of MTE is calculated wrong. I dont know why?

I calculated the values of MTE twice (another variable is masst below res1). Two different answers!!

I would be most grateful if you could help me.

Thank you

Here is my code

restart:

EPSILONE:=1:

#Digits:=10:

a[mu1]:=39.11:

b[mu1]:=533.9:

a[k1]:=7.47:

b[k1]:=0:

rhop:=5180:

rhobf:=998.2:

#gama1:=0.2;

mu1[bf]:=9.93/10000:

k1[bf]:=0.597:

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

#mu:=unapply(mu1[bf]/(1-phi(eta))^2.5,eta);

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

#mu1[bf]:=9.93/10000:

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

mu(eta)/mu1[bf];

#eq1:=diff((mu(eta)/mu1[bf])*diff(u(eta),eta),eta)+rhost(eta)-Ha^2*u(eta);

eq1:=(a[mu1]*(diff(phi(eta), eta))+2*b[mu1]*phi(eta)*(diff(phi(eta), eta)))*(diff(u(eta), eta))+(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2)*(diff(u(eta), eta, eta))+rhost(eta)-Ha^2*u(eta);

eq2:=diff((k(eta)/k1[bf])*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):

Agama1:=<<0.1>,<0.2>,<0.3>>:

ANBT:=<<0.2>,<1>,<10>>:

Aphiv:=<<0.02>,<0.06>,<0.1>>:

lambda:=0;

Ha:=0;

#NBT:=0.5:

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

#phiv:=0.02:

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

eq1;

eq2;

eq3;

res1 := dsolve(subs(NBT=1,gama1=0.3,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],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);

for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do

res := dsolve(subs(NBT=ANBT(iNBT),gama1=Agama1(igamma),phiv=Aphiv(iphi),{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],output=listprocedure,continuation=cc):

F0,F1,F2:=op(subs(subs(res),[phi(eta),u(eta),diff(T(eta),eta)])):

MTE[iNBT,iphi,igamma]:=evalf(Int((1-F0(eta)+F0(eta)*rhop/rhobf)*F1(eta), eta = 0..1)):

HTC[iNBT,iphi,igamma]:=(1+a[k1]*F0(0)+b[k1]*F0(0)^2)*F2(0):

end do;

end do;

end do;

for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do

print (convert(HTC[iNBT,iphi,igamma],float,9));

end do;

end do;

end do;

for iNBT from 1 to 3 do

for iphi from 1 to 3 do

for igamma from 1 to 3 do

print (convert(MTE[iNBT,iphi,igamma],float,9));

end do;

end do;

end do;

**Amir**