16 years, 3 days

## allvalues...

Isn't this just a case of garbage in, garbage out?
Using allvalues gives me 2 matrices:

```restart;
with(LinearAlgebra):
x := RootOf(z^2-t, z);
m := Matrix(2, (i,j) -> add(a[j, k]*((-1)^(i-1))^k*x^k, k = 0..1));
m1,m2:=allvalues(m);
L1,V1 := Eigenvectors(m1);
## and
L2,V2 := Eigenvectors(m2);```

## Lacking else...

The two outer ifs have no else, therefore if not H[i] < 2.7 or if not A[i, f] < 12 then nothing is going to happen.

So you must tell us (or yourself) what you want to do in those cases.
For clarity I wrote the code with some indentation:

```for i to n do
if H[i] < 2.7 then
if A[i, f] < 12 then
if A[i, o] < 1.2 then
Q[i, foo] := evalf(610*(A[i, o]*sqrt(H[i])*h[k]*A[i, T])^(1/2))
else
Q[i, foo] := evalf(7.8*A[i, t]+378*A[i, o]*sqrt(H[i]))
end if
end if
end if;
print(Q[i, foo])
end do;```

## Use dsolve/numeric with odeplot...

You seem to have some syntactical problems with your derivatives. Being no friend of 2D math input I made use of Kitonum's translation into the trustworthy 1D math input, known also as Maple Input.

```restart;
Digits:=15:
interface(imaginaryunit = j);
k := .5; tau := .95; mu := 0.1e-1; pi := 116.1; theta := 0.8e-2; phi := 0.25e-2; epsilon := 0.2e-2; rho := 0.5e-1; beta := 0.115e-1; chi := 0.598e-2; q := .5; eta := .2; delta := .1; alpha := 0.57e-1; p := .2; Upsilon := 1.2;
lambda := k*tau*(C(t)*Upsilon+(I)(t))/(S(t)+V(t)+C(t)+(I)(t)+R(t));
sys:={diff(C(t), t) = rho*lambda*S(t)+rho*epsilon*lambda*V(t)+(1-q)*eta*I(t)-(mu+beta+chi)*C(t), diff(I(t), t) = (1-rho)*lambda*S(t)+(1-rho)*epsilon*lambda*V(t)+chi*C(t)-(mu+alpha+eta)*I(t), diff(R(t), t) = beta*C(t)+q*eta*I(t)-(mu+delta)*R(t), diff(S(t), t) = (1-p)*pi+phi*V(t)+delta*R(t)-(mu+lambda+theta)*S(t), diff(V(t), t) = p*pi+theta*S(t)-(epsilon*lambda+mu+phi)*V(t)};
ics:= {S(0) = 8200, V(0) = 2800, C(0) = 200, (I)(0) = 210, R(0) = 200};
res:=dsolve(sys union ics, numeric);
plots:-odeplot(res,[C,I,S](t),2000..3000,numpoints=10000);
plots:-display(Array(1..5,1..1,[seq(plots:-odeplot(res,[t,F(t)],2000..3000,numpoints=10000),F=[C,I,R,S,V])]));
```

The 3 dimensional plot of C,I,S is shown here. Notice that I only plot after the initial settling has taken place. I chose 2000..3000:

Below we see that your system has 3 equilibrium points and that they all are unstable.

```RHS:=convert(evalindets(rhs~(sys),function,x->op(0,x)),list);
sol:=[solve(RHS,{C,I,R,S,V})];
pts:=map(subs,sol,[C,I,R,S,V]);
J:=unapply(VectorCalculus:-Jacobian(RHS,[C,I,R,S,V]),C,I,R,S,V);
for pp in pts do simplify(LinearAlgebra:-Eigenvalues(J(op(pp)))) end do;
```

## surd...

In maple z^(1/3) is the principal cube root of z. Thus as an example (-1)^(1/3) is equal to 1/2+I*sqrt(3)/2.
Try evalc( (-1)^(1/3) );
What you can use since you expect to get -1 from (-1)^(1/3) is surd.

Thus your plot will be produced by

```plot(surd(x^3-x^2,3),x=-3..3);
```

## Clearing some problems only...

Your code has some syntactical problems. I tried to clear them up.
Besides that the second derivatives w.r.t. eta cannot appear in the boundary conditions. I removed them.
Thus the code looks like this:

```restart;
Digits := 10;
sys_ode := 2*diff(T(eta), eta, eta, eta)+T(eta)*diff(T(eta), eta, eta) = 0;
ics := T(0) = 0, D(T)(0) = 1, D(T)(20) = 0:
sol1 := dsolve({ics, sys_ode}, numeric, output = operator);
q,w:=op(subs(sol1,[T,D(T)]));
######
PDE1 := eval([2/P * diff(g(x, eta), eta, eta)+q(eta)*diff(g(x, eta), eta)-g(x, eta)*w(eta) = 2*x*w(eta)*diff(g(x, eta), x), 2/S * diff(phi(x, eta), eta, eta)+q(eta)*diff(phi(x, eta), eta) = 2*x*w(eta)*diff(phi(x, eta), x)]);
##
subBC1 := -phi(x, 0)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
subBC2 := alpha*phi(x, 0)*sqrt(x)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
##
BC := {g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = subBC1, D[2](phi)(x, 0) = subBC2};
##
P := 1;
S := 1;
#######
pds := pdsolve(PDE1, BC, numeric, spacestep = .25);
```

Error, (in pdsolve/numeric/match_PDEs_BCs) unable to associate boundary conditions and dependent variables, system is too complex
Even when setting varepsilon=0 the same error comes up.
In an test to see if there may be other problems I tried these boundary conditions, where exp is replaced by 1:

`BCtest:={g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = -phi(x, 0), D[2](phi)(x, 0) = alpha*sqrt(x)*phi(x, 0)};`

The answer came up right away with no complaint, so it appears that the complexity of BC is the problem.

## simplify/size...

You could try simplify/size:

```length(Teller); # 4308
TS:=simplify(Teller,size);
length(TS); # 2035
whattype(TS); # `+`
nops(TS); # 5
```

## Make a procedure...

Begin by not assigning to the parameters you want to change later.

```restart;
lambda := 1.0; K__r := 1.0; S__r := .5; m := .5; M := sqrt(10.0); `&varkappa;` := .5; Omega := sqrt(5.0);
## Have omitted Gm, Gr, P__r, S__c
## Now define PDE and IBC as you did
##  ...
## Define SOL as a procedure:
SOL:=(gm,gr,pr,sc)->pdsolve(eval(PDE,{Gm=gm,Gr=gr,P__r=pr,S__c=sc}), IBC, numeric, spacestep = 0.1e-1);
## Test:
SOL(5,6,0.71,0.22):-plot(t = .3, color = red);
SOL(5,6,0.71,0.1):-plot(t = 1,color = red);
```

Adding this code you can get an animation;

```PL1:=sc->SOL(5,6,0.71,sc):-plot(t = 1,color = red);
animate(PL1,[sc],sc=0.01..1);
```

## Syntax...

I'm assuming that you meant sys:= {...} and not sys*{...}.
Thus we have:

```sys:={A2+D2 = 0, B1*sin(192*K1)+D1 = 0, 3.383*B2*K1+C2 = 0, B1*K1^2*sin(192*K1)+11.444689*A2*K1^2*cos(568.344*K1)+11.444689*B2*K1^2*sin(568.344*K1) = 0, A2*cos(568.344*K1)+B2*sin(568.344*K1)+168*C2+D2 = 0, -3.383*A2*K1*sin(568.344*K1)+3.383*B2*K1*cos(568.344*K1)+C2-B1*K1*cos(192*K1) = 0};
##You cannot use fsolve if you want to keep K1 unassigned.
solve(sys,{A2, B1, B2, C2, D1, D2});
```

You just get the trivial solution {A2 = 0., B1 = 0., B2 = 0., C2 = 0., D1 = 0., D2 = 0.}.
## You may ask: For which values of K1 are there nontrivial solutions?
Notice that your equations are linear in A2, B1, B2, C2, D1, D2.
To look into that use a little linear algebra:

```A,zz:=LinearAlgebra:-GenerateMatrix(sys,[A2, B1, B2, C2, D1, D2]);
dt:=LinearAlgebra:-Determinant(A); #Must be zero to get nontrivial solutions
plot(dt,K1=-.1..0.1);
plot(dt,K1=-.01..0.01);
```

We must expect infinitely many values of K1 with that property.
# Since dt-eval(dt,K1=-K1) evaluates to zero, you need only look at K1 >= 0.
To find the corresponding solutions for A2, B1, B2, C2, D1, D2 you could do:

```K1a:=fsolve(dt=0,K1=0.006); #Example only
G:=LinearAlgebra:-GaussianElimination(eval(A,K1=K1a));
G0:=fnormal~(G); #When dt=0 the rank < full. Removing roundoff.
LinearAlgebra:-LinearSolve(G0,zz,free=t);
```

## convert/exp...

```S:=cos(440*2*Pi*t)+cos(554*2*Pi*t)+cos(659*2*Pi*t);
convert(S,exp);
```

I changed pi to Pi.

## Right clicking...

This probably won't be satisfactory, but I tried without using plotsetup. I right clicked on the plot and exported as a jpg file:

```plots:-logplot(x^2-3, x = 0 .. 100, font = [bold, "TimesNewRoman", 30],size=[1500,1100]);
```

The plot:

## _Fn...

The functions named _F1, _F2, ... are meant as arbitrary functions. Thus an expression like _F1(t+x) means any expression depending on the sum t+x only. Nothing more is meant (except that appropriate smoothness requirements are implicitly assumed).

## A few changes...

You need initial conditions for the derivative too. Your ode is called ode1, not as you have ode.
To get arrows you must rewrite as a system of first order:

```restart;
ode1 := diff(x(t), t, t) = (5*9.80665)*sin((1/6)*Pi)-(10*(10-sqrt(x(t)^2+25)))*x(t)/sqrt(x(t)^2+25)-(diff(x(t), t));
ics:=[seq([x(0) = (1/4)*k,D(x)(0)=k/2], k = -20 .. 20)]; # I chose D(x)(0)=k/2
DEtools[DEplot](ode1, x(t), t = -2 .. 2, ics, x=-7..20, stepsize = 0.5e-1, linecolour = red);
ode2:=diff(x(t),t)=y(t);
sys := {ode2,subs(ode2,ode1)};
ics_xy:=subs(D(x)=y,ics);
DEtools[DEplot](sys, [x(t),y(t)], t = -2 .. 2, ics_xy, x=-7..15,y=-7..20,color = blue, stepsize = 0.5e-1, linecolour = red,arrows=medium);
```

The last plot:

## Cleaning up the solution...

As we all agree this is just a "cleaning up" business. There is nothing wrong with the solution.
I just tried to make a procedure CleanUp. It assumes that the input is a solution to a linear ode and comes from dsolve.
It hasn't been tested thoroughly.
For what it is worth here it is:

```CleanUp:=proc(sol::equation) local t,cs,Hom,Inhom,S,tm;
if not lhs(sol)::function(name) then error "This is not an output from dsolve" end if;
t:=op(1,lhs(sol));
cs:=indets(sol,'suffixed'(_C,posint));
if cs={} then return sol end if;
Hom,Inhom:=selectremove(has,rhs(sol),cs);
if not type(Inhom,`+`) then return sol end if;
S:={};
for tm in Inhom do
if solve('identity'(Hom=tm,t))<>NULL then S:=S union {tm} end if
end do;
Inhom:=subs(S=~0,Inhom);
lhs(sol)=Hom+Inhom
end proc:
```

## You found it...

Try this:

```Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, method=newtoncotes[open,4], output=information);
```

There is a discrepancy, however, between what is written in the help page
under Options and what is actually the case.
Whether to apply the adaptive version of the specified method. The following methods are available when adaptive is set to true:
Newton-Cotes Rule, open or closed, with a user-specified order"

but if you try

```Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, adaptive=true,method=newtoncotes[open,4], output=information);
```

then you will get the error message

Error, (in Student:-NumericalAnalysis:-Quadrature) the newtoncotes[open, 4] method is not available when the adaptive option is specified

##
Note: I have submitted an SCR.

## Solving numerically...

I changed pi into Pi and made an assignment to K in your worksheet.
I chose as an example to make all parameters in the system equal to 1 and all initial values equal to zero.
Change that as appropriate.
MaplePrimes18-03-07odesys.mw
To get the plots of all 4 separately you can use this code:

```PLTS:=[seq(plots:-odeplot(res,[t,F(t)],0..7,numpoints=10000),F=[V1,V2,V3,V4])]:
plots:-display(Array(1..4,1..1,PLTS));
```

The plot at the end showing V1 and V2 only (and together):

 First 12 13 14 15 16 17 18 Last Page 14 of 149
﻿