Following previous question at

http://www.mapleprimes.com/questions/149581-Improve-Algorithm-Dsolve

and also

http://www.mapleprimes.com/questions/149243-BVP-With-Constraining-Integrals

I wrote the following code

***********************

restart:

gama1:=0:

phi0:=0.00789:

rhocu:=2/(1-zet^2)*int((1-eta)*rho(eta)*c(eta)*u(eta),eta=0..1-zet):

eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])+((1/(eta-1)+1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta)):

eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*(2/(1-zet^2)*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)-k(eta)/k1[w]/(1-eta)*diff(T(eta),eta) )):

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

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:=0.5:

#phi(0):=1:

#u(0):=0:

phi1[w]:=phi0:

N[bt]:=0.2:

mu1[w]:=mu(0):

k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,u(0)=0,eq1):

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

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

p:=proc(pp2) global res,F0,F1,F2:

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

res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure):

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

evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2*10000:

end proc:

s1:=Student:-NumericalAnalysis:-Secant(p(pp2),pp2=[6,7],tolerance=1e-6);

HFloat(6.600456858832996)

p2:=%:

ruu:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet))):

phb:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet))) :

TTb:=evalf(2/(1-zet^2)*(Int((1-eta)*F2(eta),eta=0..1-zet))):

rhouu:=evalf(2/(1-zet^2)*(Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))):

with(plots):

res(parameters=[R0,R1]):

odeplot(res,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..zet);

*************************************

as you can see at the second line of the code, the value of phi0:=0.00789. however, I want to modify the code in a way that phi0 is calculated with the following addition constraint

evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)))-0.02=0

I would be most grateful if you could help me in this problem.

Thanks for your attention in advance

**Amir**