16 years, 0 days

## Functions...

It seems that algsubs won't touch attempts to replace functions like int, Int, sin, or f in f(x):

```restart;
expr:=int(f(x),x);
algsubs(f(x)=g(x),expr); # Works
algsubs(f=g,expr); # Doesn't do anything
algsubs(sin=cos,sin(x)); # Doesn't do anything
algsubs(x=y,expr); # Works
```

## Infinitely many solutions...

That problem with y(0) = 0 has infinitely many solutions (a well known fact):

```restart;
ode := diff(y(x),x) = 2*sqrt(y(x));
p:=seq(piecewise(x<n,0,(x-n)^2),n=0..4);
plot([p],x=-1..5);
seq(odetest(y(x)=q,ode),q=p);
```

## Introduce extra equation(s)...

With 17 equations we may assume that your solution is numerical.

I give here a very simple example: Just one ode and it can be solved symbolically as well as numerically:

```restart;
ode1:=diff(x(t),t)=-x(t)+sin(t);
ic1:=x(0)=1;
res1:=dsolve({ode1,ic1},numeric); # The numerical solution for x
plots:-odeplot(res1,[t,x(t)],0..20);
### Now introduce the integral of x which is called y:
INT:=y(t)=int(x(s),s=0..t);
### Differentiate that to get an extra ode:
ode2:=diff(INT,t);
ic2:=eval(INT,t=0); # The initial value for y
res2:=dsolve({ode1,ode2,ic1,ic2},numeric); # Numerical solutions for x and y
plots:-odeplot(res2,[[t,x(t)],[t,y(t)]],0..20);
### Now for comparison we solve exactly (symbolically):
sol1:=dsolve({ode1,ic1});
map(int,sol1,t=0..T);
sol2:=y(t)=subs(T=t,rhs(%));
plot([rhs(sol1),rhs(sol2)],t=0..20);
### The difference between the numerical and the exact y(t):
plots:-odeplot(res2,[t,y(t)-rhs(sol2)],0..20);
```

PS. I see that Carl Love beat me to it, but I shall leave my version. It doesn't hurt with several examples especially when you are new to Maple.

## Events?...

Could you not use events?
Here is a suggestion:

```sysE := {
diff(x(t), t) = v(t),
diff(v(t), t) = t - (b(t)*f+(1-b(t))*(-f)),
x(0) = 0,
v(0) = 0,
b(0) = 0
};
solE := dsolve(sysE, numeric, parameters=[f],discrete_variables=[b(t)::boolean], events=[[v(t)=0,b(t)=1-b(t)]]);
solE(parameters=[1]):
solE(1);
odeplot(solE,[t,v(t)],0..100);
```

## Use assumptions...

In the revised version of U it pays to reverse the order of integration:

```restart;
#HAM := [1, 2, 3, 4, 5, 6];
U := q^6*sin(p);
alpha := 1/3;
## Notice that the order of integration has been reversed:
for i from 1 to 6 do
rho := i;
f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), p = 0 .. x), q = 0 .. y) assuming x>0,y>0
end do;
```

## A result (finally)...

 > restart;
 > sys:=P*diff(f(eta), eta, eta, eta)/Q + lambda*theta(eta)*(S0 + S1*Qc*theta(eta)) + (n + 1)/2*(f(eta) + g(eta))*diff(f(eta), eta, eta) - n*diff(f(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, P*diff(g(eta), eta, eta, eta)/Q + (n + 1)/2*(f(eta) + g(eta))*diff(g(eta), eta, eta) - n*diff(g(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, (S2 - 2*R)*diff(theta(eta), eta, eta)/(S4*Pr) + (3*R)/(2*S4*Pr*(thetaw - 1))*diff((1 + (thetaw - 1)*theta(eta))^2, eta, eta) + (n + 1)/2*(f(eta) + g(eta))*diff(theta(eta), eta) = 0, f(0) = 0, g(0) = 0, theta(0) = 1, D(f)(0) = 1, D(g)(0) = c, D(f)(3) = 0, D(g)(3) = 0, theta(3) = 0;
 > sys[3]; # Term 2 has a removable singularity at thetaw=1
 > sys3:=map(normal,lhs(sys[3]))=0; # No singularity at thetaw=1
 > isolate(sys3,diff(theta(eta),eta,eta));

When R > 0, theta(eta) > 0, and thetaw > 1 we see that the denominator is positive (in fact grater than R+S2).

 > SYS:=sys[1..2],sys3,sys[4..-1];
 > P := 0.7952865489: Q := 1.239340734: S4 := 1.007143462:S0 := 0.8330195568:S1 := 0.8330195568: S2 := 1.087762421: Pr := 0.0157*2415/0.252: n:=3: lambda:=20: Qc:=2: c:=0.5:
 > for s in SYS do s end do; # Viewing the system
 > Digits:=15:

Boundary value problems of this kind are notoriously difficult, so we begin with simply providing simple values for thetaw and R:

 > res:=dsolve(eval({SYS},{thetaw=1.5,R=0.2}),numeric,initmesh=256,maxmesh=1024);
 > res(0);
 > dth0:=subs(res(0),diff(theta(eta),eta));
 > plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3);

This procedure will be used for animation:

 > resplot:=proc(thetaw0,R0) local res;    res:=dsolve(eval({SYS},{thetaw=thetaw0,R=R0}),numeric,_rest);    plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3,legend=[f,g,theta]) end proc;
 > resplot(1.5,0.2,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]);
 > plots:-animate(resplot,[1.5,R,initmesh=512,maxmesh=2048,        approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],R=0.1..2,trace=24);
 > plots:-animate(resplot,[thetaw,0.2,initmesh=512,maxmesh=2048,               approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],thetaw=1..10,trace=24);

We notice that the general shape of the solutions doesn't change much with R and thetaw.

The following procedure takes values of thetaw and D(theta)(0) as input and outputs the corresponding value for R.

Notice the extra boundary condition D(theta)(0)=Dth0:

 > Rproc:=proc(thetaw0,Dth0) local res;    if not [thetaw0,Dth0]::list(realcons) then return 'procname(_passed)' end if;    res:=dsolve(eval({SYS,D(theta)(0)=Dth0},thetaw=thetaw0),numeric,_rest);    subs(res(0),R) end proc;

Here we should expect to see the result 0.2 for R even though the R value given in approxsoln is R = 1, and indeed we do:

 > Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);
 > Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=2]);
 > Rproc(1.5,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);
 > plot(      Rproc(thetaw,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),      thetaw=1..2,labels=[thetaw,R]);
 >

Here as we saw we have D(theta)(0) = -1 (the second argument to Rproc):

This code:

```plot3d(Rproc(thetaw,Dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
thetaw=1..2,Dth0=-3..-1,labels=[thetaw,Dth0,R]);
```

produces (after a while):