Dear Colleges

I have a problem with the following code. As you can see, procedure Q1 converges but I couldn't get the resutls from Q2.

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

Sincerely yours

Amir

restart;

Eq1:=diff(f(x),x$3)+diff(f(x),x$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x$2)^2*x^2,x$1);

Eq2:=diff(g(x),x$3)+diff(g(x),x$2)*g(x)+c*a^2*sqrt(2*reynolds)*diff(diff(g(x),x$2)^2*x,x$1);

eq1:=isolate(Eq1,diff(f(x),x,x,x));

eq2:=subs(g=f,isolate(Eq2,diff(g(x),x,x,x)));

EQ:=diff(f(x),x,x,x)=piecewise(x<c*0.1,rhs(eq1),rhs(eq2));

Eq11:=diff(theta(x),x$2)+pr*diff(theta(x),x$1)*f(x)+pr/prt*b^2*sqrt(2*reynolds)*diff(diff(f(x),x$2)*diff(theta(x),x$1)*x^2,x$1);

Eq22:=diff(g(x),x$2)+pr*diff(g(x),x$1)*f(x)+pr/prt*a^2*c*sqrt(2*reynolds)*diff(diff(f(x),x$2)*diff(g(x),x$1)*x^1,x$1);

eq11:=isolate(Eq11,diff(theta(x),x,x));

eq22:=subs(g=theta,isolate(Eq22,diff(g(x),x,x)));

EQT:=diff(theta(x),x,x)=piecewise(x<c*0.1,rhs(eq11),rhs(eq22));

EQT1a:=eval(EQT,EQ):

EQT2:=eval(EQT1a,{f(x)=G0(x),diff(f(x),x)=G1(x),diff(f(x),x,x)=G2(x)}):

bd:=c;

a:=0.13:

b:=0.41:

pr:=1;

prt:=0.86;

reynolds:=12734151.135786774055543653356602; #10^6; #1.125*10^8:

c:=88.419896050808975395120916434619:

;

Q:=proc(pp2) local res,F0,F1,F2;

print(pp2);

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

res:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=pp2},numeric,output=listprocedure);

F0,F1,F2:=op(subs(subs(res),[f(x),diff(f(x),x),diff(f(x),x,x)])):

F1(bd)-1;

end proc;

fsolve(Q(pp2)=0,pp2=(0..1002));

se:=%;

res2:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=se},numeric,output=listprocedure):

G0,G1,G2:=op(subs(subs(res2),[f(x),diff(f(x),x),diff(f(x),x,x)])):

plots:-odeplot(res2,[seq([x,diff(f(x),[x$i])],i=1..1)],0..c);

Q2:=proc(rr2) local solT,T0,T1;

print(rr2);

if not type(rr2,numeric) then return 'procname(_passed)' end if:

solT:=dsolve({EQT2,theta(0)=1,D(theta)(0)=-rr2},numeric,known=[G0,G1,G2],output=listprocedure):

T0,T1:=op(subs(subs(res),[theta(x),diff(theta(x),x)])):

T0(bd);

end proc;

fsolve(Q2(rr2)=0,rr2=(0..100));

shib:=%;

sol:=dsolve({EQT2,theta(0)=1,D(theta)(0)=-shib},numeric,known=[G0,G1,G2],output=listprocedure):

plots:-odeplot(sol,[x,theta(x)],0..c);

#fsolve(Q2(pp3)=0,pp3=-2..2):

**Amir**