Hi,

according to my previous question

http://www.mapleprimes.com/questions/201435-ODE-With-Constraint

I wrote the following code. at first, the code solve the equation for f and when it slves that, I want to solve Theta in such a way that use the values of f in previous calculation. I use the command 'known' but i couldnt find thesolution

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

Thanks for your attentions in advance

restart; # Notice that Restart (capital R) has no effect (to catch that use semicolon, not colon)

a:=0.13:

b:=0.41:

reynolds:=1.125*10^7;

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

c:=75:

;

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(c)-1;

end proc;

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

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=0..2)],0..2); #This plots from and past 0.1*c

pr:=1;

prt:=0.89;

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(f(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));

QT:=proc(pp3) local res3,theta0,theta1;

print(pp3);

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

res3:=dsolve({EQT,theta(0)=1,D(theta)(0)=pp3,known=f},numeric,output=listprocedure);

theta0,theta1:=op(subs(subs(res),[theta(x),diff(theta(x),x)])):

theta0(c);

end proc;

fsolve(QT(pp3)=0,pp3=(0..200));

res3(0);

**Amir**