## 11764 Reputation

15 years, 325 days

## taylor...

Maybe you want the coefficients in a taylor expansion in p?
I'm using mtaylor instead of taylor for convenience since it omits the big O term:
So you can do:

```T_HPMEq1:=mtaylor(HPMEq1,p,N+1):
T_HPMEq2:=mtaylor(HPMEq2,p,N+1):
T_HPMEq3:=mtaylor(HPMEq3,p,N+1):
#
for i from 0 to N do
equ1[i] := coeff(T_HPMEq1, p, i) = 0;
equ2[i] := coeff(T_HPMEq2, p, i) = 0;
equ3[i] := coeff(T_HPMEq3, p, i) = 0
end do;
```

## DEtools[remove_RootOf]...

Since odetest handles the implicit given solution nicely and regardless of the name or value of the arbitrary constant (see my reply to Carl Love) another workaround is to convert the "explicit" solution to an implicit one. I use quotes since a solution having a RootOf can hardly be called a genuinely explicit one.

```restart;
ode := diff(y(x), x) = x*ln(y(x)):
sol:=dsolve(ode, implicit);
odetest(sol,ode); #OK
odetest(subs(_C1=C1,sol),ode);#OK
odetest(subs(_C1=1,sol),ode); #OK
explicit_sol:=op(solve(sol,{y(x)}));
DEtools[remove_RootOf](subs(_C1=C1,explicit_sol));
odetest(%,ode);
DEtools[remove_RootOf](subs(_C1=1,explicit_sol));
odetest(%,ode);
```

## cat...

Use cat.

```x1, x2, x3, x4:=1,2,3,4;
for i from 1 to 4 do unassign(cat(x,i)) end do;
x1,x2,x3,x4;
```

As for || see the help page ?|| where it says:
"Note: The use of this operator is discouraged. Where possible the function cat should be used."

## _EnvDSNumericSaveDigits...

I shall just address the warning about _EnvDSNumericSaveDigits. That is curious and I haven't seen that before.
It seems to come from `dsolve/numeric` or rather from one ot its auxiliary procedures.
Notice that the warning disappears if you make an explicit integer assignment to _EnvDSNumericSaveDigits such as

`_EnvDSNumericSaveDigits:=15:`

Since the environmental variable _EnvDSNumericSaveDigits is meant to save the old value of Digits any positive integer will make the warning go away.
The warning must come up in in your worksheet because somehow _EnvDSNumericSaveDigits avoided getting assigned to in your worksheet. The warning doesn't come up by a simple call to any of the two procedures, only in the fsolve context (or so it seems).
### Furthermore if you remove the range argument then the warning disappears.
Here is a much simpler example where the same warning comes up:

```restart;
ode := diff(x(t), t) = x(t);
t1:=2:
h := proc (x0)
local X;
X := eval(x(t), dsolve({ode, x(0) = x0}, numeric,range=0..t1, output = listprocedure));
X(t1)-1
end proc:
h(0);
fsolve(h);
h(%);
```

As you see the result is correct, but the warning is confusing.
Try removing the range argument and the warning is gone.
I would have done this example using the parameters option in this way:

```restart;
ode := diff(x(t), t) = x(t);
t1:=2:
res:=dsolve({ode, x(0) = x0}, numeric,range=0..t1, output = listprocedure,parameters=[x0]);
X:=subs(res,x(t));
h := proc (x0)
res(parameters=[x0]);
X(t1)-1
end proc:
h(0);
fsolve(h);
h(%);
```

No warnings. You can remove the range argument (and still no warnings).
I shall submit an SCR (bug report) about the warning as it is confusing to the user.

## Save both...

If you save both as .mpl files and use the complete path and then use a read statement to read main_module.mpl your setup with \$include should work.
I just tried:

```main_module:=module()
option package;
local C;
local sub_module;
export main_entry;

C:=99; #see if this can be "seen" from child module

\$include "F:/MapleDiverse/MaplePrimes/sub_module.mpl"

main_entry :=proc()
print("in main_module:-main_entry()"):
sub_module:-main_entry();
end proc;
end module;```

and

```#sub_module.mpl
sub_module:= module()
export main_entry;
main_entry :=proc()
print("In main_module:-sub_module:-main_entry(), C=",C):
end proc;
end module;```

The files are in my case both in the folder "F:/MapleDiverse/MaplePrimes".
Now in Maple use the read statement

```restart;
## Test:
main_module:-main_entry();```

## A symbolic solution...

It appears that Maple 2018 can find a symbolic solution. It is big, but integrating the first ode w.r.t. x and using ics helps and is done here:
Note added: The code below assumes that Q is a constant. If Q depends on x the method may still work. If not use numerical solution as Carl does.

```restart;
sys_ode:=    2*C__5*diff(w(x),x) - C__1*diff(u(x),x\$2) - 2*C__4*diff(w(x),x\$3) = Q,
2*C__2*u(x) - C__1*diff(w(x),x\$3) - 2*C__3*diff(u(x),x\$2) = 0;
ics:= w(0) = 0, D(u)(0) = 1, (D@@2)(w)(0) = 0, D(u)(100) = 1, (D@@2)(w)(100) = 2;
###
map(int,sys_ode[1],x=0..x,continuous);
ode1:=eval(%,{ics});
ode2:=sys_ode[2];
### Now solve symbolically:
res:=dsolve({ode1,ode2,ics}):
U,W:=op(subs(res,[u(x),w(x)])):
Digits:=200:
## Example plot:
plot(eval(U,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);

```

The plot of W is done similarly:

```plot(eval(W,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);
```

## Formal parameter...

In procB you have the formal parameter the_equation. That appears in the procedure body of procB in the print statement as well as in the call to procA. Wen procB is called with the equation y=3 the formal parameter the_equation will be replaced by y=3 in all its occurrences in the body of procB. Thus the call to procA will be procA('y=3' = (y=3)), so that the default NULL will be used by procA.

`procA(':-the_equation'=the_equation);`

## Inevitable...

Your sys1 will necessarily have a singularity for some t1 < tini.
The reason is that as long as a(t) > 0 it follows that a(t) is increasing. Going back in time a(t) will decrease to 0 at some finite value of t (called t1 above). Here a'(t1) = infinity.
To prove that that is what happens is rather straightforward. Your ode is of the form

`ode:=diff(a(t),t) = k*a(t)*sqrt(1/a(t)^3 + m); # k>0, m>0.`

with a(t0) = a0 > 0.

As long as a(t) > 0 we have

`diff(a(t),t) > k*a(t)*sqrt(1/a(t)^3) = k/sqrt(a(t))`

Rewrite this to get

`sqrt(a(t))*diff(a(t),t) > k`

Integrate from t to t0 (with t < t0) to get

`2/3*(a0^(3/2) - a(t)^(3/2)) > k*(t0-t)`

from which follows that

`t > t0 - a0^(3/2)/k*(2/3)`

since we are assuming that a(t) > 0. Thus a(t) will necessarily approach 0 from above at a finite value of t < t0.

## expand...

Try this where I have replaced the letter i by the letter I:

```restart;
A:=Matrix(6,(i,j)->cos(2*Pi/6*(i-1)*(j-1))-I*sin(2*Pi/6*(i-1)*(j-1)));
B:=conjugate~(A);
expand~(A.B);
```

## Eliminate one of the two...

If you eliminate one of CNO and CNH3 you can easily solve the system. From your two odes it follows that the difference between CNO and CNH3 is constant, thus in the case presented equal to zero.

```restart;
r := 2900*exp(-53300/(R*T))*CNO(t)^0.62*CNH3(t)^(-0.05);
ode := diff(CNO(t), t) = -r: #(negative because they are decaying with time)
ode2 := diff(CNH3(t), t) = -r:
ics := CNO(1020) = 1.6, CNH3(1020) = 1.6; #(sets up known initial conditions)
ode := subs(R = 8.314, T = 473, ode):
ode2 := subs(R = 8.314, T = 473, ode2):
#sys_ode := (ode, ode2) ;
odeS:=subs(CNH3(t)=CNO(t),ode);
#dsolve({odeS,CNO(1020)=1.6}); #??? Why does it take so long (or forever?)
sol:=dsolve(odeS);
eval[recurse](sol,{t=1020,CNO(1020)=1.6});
solve(%,{_C1});
eval(sol,%);
res:=solve(%,CNO(t)); # The solution for CNO(t)
t0:=solve(res=0,t);
plot(res,t=1020..t0);
```

The solution becomes 0 in finite time t0 (approximately 1775.4) due to the power of the right hand side of odeS being positive and less than 1. After that the solution continues as zero.

If the initial conditions are different it appears that you run in to an integral that Maple cannot find. A numerical solution may then be used:

```ics2 := CNO(1020) = 1.6, CNH3(1020) = 1; #Example of different initial conditions
# Now CNH3(t)-CN0(t) = constant = 1 - 1.6 = -0.6
odeS2:=subs(CNH3(t)=CNO(t)-0.6,ode);
res:=dsolve({odeS,CNO(1020)=1.6},numeric);
plots:-odeplot(res,[t,CNO(t)],1020..2000);
```

Notice though that since CNH3(t) = CNO(t) - 0.6 we must stop the integration when CNH3(t) = 0.
You could just plot CNH3(t) like this:

plots:-odeplot(res,[t,CNO(t)-0.6],1020..2000);

An alternative:

```resE:=dsolve({odeS,CNO(1020)=1.6},numeric,events=[[CNO(t)-0.6,halt]]);
resE(2000);```

You will be told that integration stopped at t = 1279.9381.

## Color...

Make the following addition and change:

```col:=["black", "blue", "red", "pink"];
for k to 4 do
R := dsolve(eval({bc, sys}, Ha = L[k]), fcns, type = numeric, AP);
AP := approxsoln = R;
p1u[k] := odeplot(R, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], style = line, color = col[k])
end do:
```

## A way, but with problems...

Begin by revising the definition of F and then use a revised version of your second method.
## I revised my answer too  :-)
## Introduce Y'' as Yd2(x) (the easiest solution in this case, I think):

```F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
Y1 := subs(F,Y(x));
Y2 := subs(F,diff(Y(x), x));
Y3 := subs(F,Yd2(x));

```

The remaining ones ( R(x), mu(x) , and mu'(x) ) are found quite similar to Y(x), and Y'(x).

## There appears to be singularity problems.

```F := dsolve({cond, sys}, numeric, output = listprocedure); #No Yd2
plots:-odeplot(F,[x,Y(x)],0..Pi);
```

You seem to be interested in x = 0..Pi, but you get from odeplot:

Warning, cannot evaluate the solution further left of 2.6226594, probably a singularity

Furthermore when trying to include diff(Y(x),x,x) = Yd2(x):

```F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
plots:-odeplot(F,[x,Y(x)],0..Pi);
```

then we get an error I don't remember having ever seen before:

Error, (in dsolve/numeric/SC/IVPrun) invalid reversal of integration direction

Going the other way is OK at least until 3.66:

```plots:-odeplot(F,[x,Y(x)],Pi..3.66);
```

# When I use the second version of F (i.e. the one including Yd2) and try F(3) then I get this error:

Error, (in unknown) cannot evaluate the solution further left of 3.6605260, probably a singularity

which is surprising to me because who told F to go from Pi to 3.6605260 before going to 3?

## I notice that the weird error disappears if cot(x) is replaced by 1. Singularity problems don't disappear though.
###############
The problem (unrelated to other singularities) appears to have to do with the initial singularity at x = Pi where cot(x) is infinite (undefined):
Try

```F := dsolve({cond, sys}, numeric, output = listprocedure,abserr=1e-15,relerr=1e-13);
F(3);
```

and you get the error message:

Error, (in unknown) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

## Probably indeed a singularity!...

Actually, in Maple 2018 I get the complaint:

Warning, cannot evaluate the solution further right of 3.4774930, maxfun limit exceeded (see ?dsolve,maxfun for details)
So since maxfun can be set, I did this:

`res := dsolve(sys union ics, numeric, maxfun=10^6);`

But the complaint just changed the t-value slightly to 3.4775188.
There may be an actual singularity though.
Try this:

```sys:=solve(odes,diff~({x,y,z}(t),t)); # Solve for the derivatives
## Then do
plots:-odeplot(res, [t,rhs(sys[1])], Pi .. 4,caption="The derivative of x(t)");

```

What I see is this:

Using method=rosenbrock gives me the error you report:

`res := dsolve(sys union ics, numeric, maxfun=10^6,method=rosenbrock);`

## Function definition...

The correct way to define a function is like this:

```T:=x->1/2*x*sqrt(-x^2+(60-x)^2);
##
solve(diff(T(x),x)=0,x);
##
T(20);
```

## Calculus1:-ApproximateInt, NumericalAnal...

If you haven't already done so, you should take a look at Student:-Calculus1:-ApproximateInt and Student:-NumericalAnalysis:-Quadrature at least for comparison with your own code.

```restart;
fx:=x^2*sqrt(25-x^2);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=trapezoid,partition=100);
evalf(%);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=midpoint,partition=100);
evalf(%);
```Student:-NumericalAnalysis:-Quadrature(fx,x=0..5,method=newtoncotes[open,0], partition=100);