Question: Wrong calculations of the results

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