A:=Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]):
B := Matrix(ListTools:-Reverse(convert(A, listlist)));

I understood the question that several separately created animations should simultaneously (in parallel) be reproduced on the same plot. To obtain the desired ranges in both axes, use the  view  option. Example below:

A:=plots:-animate(plot,[x^2, x=-2..a, color=red], a=-2..2, frames=60):
B:=plots:-animate(plot,[sin(2*x), x=0..a, color=blue], a=0..2*Pi, frames=60):
plots:-display(A, B, view=[-2..2*Pi, -1..4]);

If you are interested in various ways of combining several animations (parallel and sequential), then look at this post  https://www.mapleprimes.com/posts/207840-Combinations-Of-Multiple-Animations

Maple does not support the solution of symbolic matrix equations, even in the simplest case  A.X=B  when A is a invertible matrix. Similar equations can be solved in  LinearAlgebra  package only for specific  matrices  A  and  B .

restart;
N:=5:
A:=Matrix(N):
for i from 1 to N do
for j from 1 to N do
A[i,j]:=1/(i*j);
od: od:
A;

Of course, the same can be done in short:

A:=Matrix(5, (i,j)->1/(i*j));

restart;
simplify([a^4 + b^4 + c^4, a^5 + b^5 + c^5], {a+b+c=1, a^2+b^2+c^2=2, a^3+b^3+c^3=3});
solve({a+b+c=1, a^2+b^2+c^2=2, a^3+b^3+c^3=3});
evalf(%);

Here are the solutions:

# Question 1:
f:=x->x^3-2*x;
a:=0:  b:=2:
solve(f(b)-f(a)=D(f)(c)*(b-a), c);
evalf(%);

# Question 2:
g:=x->surd((x-3)^2, 3);
a:=-3:  b:=4:
solve(g(b)-g(a)=D(g)(c)*(b-a), c);
evalf(%);

In the first example of the two solutions obtained, we should select only the positive root.

Note that in the second example the point  c  does not belong to the interval   [-3,4] . The reason is that the function  g  is not differentiable at the point  x=3 . See the plots:

plot(f, 0..2);
plot(g, -3..4);

An arbitrary point inside the triangle you can get as a convex linear combination of its vertices. See  https://en.wikipedia.org/wiki/Convex_combination

For example:

restart; 

A := [1, 2, 3]; 
DD := [-2, 1, 0]; 
T := [1/2, 1/2, 3/2];

P:=1/3*A+1/3*DD+1/3*T;

P is the intersection point of the medians of this triangle.

In the real domain, this identity is also true for  x< 0 . Since  coth(x)*tanh(x)=1 , we will simply replace  1  in the expression  f  with  coth(x)^(1/3)*tanh(x)^(1/3) :

restart;
	
f:=(coth(x)^(1/3)-tanh(x)^(1/3))*(coth(x)^(2/3)+tanh(x)^(2/3)+coth(x)^(1/3)*tanh(x)^(1/3));

expand(f);
simplify(%);

restart;
Expr:=A*exp(-m*x)*cos(y)+B*exp(-m*x)*sin(y):
solve(identity(Expr,y), {A,B});
                                          {A = 0, B = 0}

The shape of the container can be very different. We restrict ourselves to the containers in the shape of a cuboid, and we will assume that length>=width>=height. The task is reduced to the factorization of the number 320 by 3 factors taking into account the conditions above. For this we use the procedure  Factoring  from this post  https://www.mapleprimes.com/posts/141668-Partitions-Of-A-Natural-Number-Into-Factors

ListTools:-Reverse~([map(t->[1,t[]],Factoring(320,2))[],Factoring(320,3)[]]);

[[160, 2, 1], [80, 4, 1], [64, 5, 1], [40, 8, 1], [32, 10, 1], [20, 16, 1], [80, 2, 2], [40, 4, 2], [32, 5, 2], [20, 8, 2], [16, 10, 2], [20, 4, 4], [16, 5, 4], [10, 8, 4], [8, 8, 5]]

This choise  [320,1,1]  is omitted for aesthetic reasons. Perhaps the first few variants (of those that remain) should also be discarded.

There are some problems with this code.

Here is an example with a hemisphere. In the case 1, there are many extraneous lines, in the case 2 the code does not work at all:

L:=[cos(x), sin(x)*cos(y), sin(x)*sin(y)]:
plot3d(L,x=0..Pi/2, y=0..2*Pi, scaling=constrained); # The surface 
ContoursWithLabels(L, 0..Pi/2, 0..2*Pi, {0.3,0.7,0.9}, [scaling=constrained]); # The case 1
ContoursWithLabels(L, 0..Pi/2, 0..2*Pi, [scaling=constrained]); # The case 2
 

You can get a good approximation period using best approximations  sqrt(2)  with rational numbers. See  https://en.wikipedia.org/wiki/Square_root_of_2  For example using  sqrt(2)=577/408  we get good approximation period  T=204*Pi .

Compare (in the third plot, the previous 2 ones are superimposed on each other.):

A:=plot((100+100*cos(6*t)+200 *cos(12*sqrt(2)*t))^2, t=0..3, color=red, size=[800,200]);  
B:=plot((100+100*cos(6*t)+200 *cos(12*sqrt(2)*t))^2, t=204*Pi..204*Pi+3, color=blue,size=[800,200]);
plots:-display(A,plottools:-translate(B,-204*Pi,0));  

There is another simple way to prevent the immediate calculation of a trigonometric function. If the argument contains  Pi , then just replace  Pi  with  pi . On the output  pi  will look exactly the same as  Pi . When using the  % symbol with the name of a constant or function, they will be displayed in gray, which is not always desirable:

A:=%cos(2*Pi/3)+I*%sin(2*Pi/3);
B:=cos(2*pi/3)+I*sin(2*pi/3);
eval(B, pi=Pi);
                              

Do not use the name and the indexed name with the same symbol in the same document. In your code these are  c  and  c[1], c[2], c[3] . Take a look at what might happen in this case:

restart;
c:=6:
c[1];
                              

I replaced  с  with  С , but then a new error appears. Probably the model itself is incorrect.

restart;
with(plots):
Q := 1000: a := .9: b := 0.8e-3: mu := 0.247e-2: k := .2: y := 0.2e-1: e := .5: g := .5: T := 8: n := 100:
sigma[1] := 0.5e-1: sigma[2] := 0.9e-1: alpha[1] := 0.52e-2: alpha[2] := 0.52e-2: delta[1] := .8: delta[2] := .904: delta[3] := .8: c[1] := 50: c[2] := 250: c[3] := 50: w[1] := 140: w[2] := 130: w[3] := 150: w[4] := 160:
u[1] := min(max(0, z), 1): z := (a*i[p](t)*(i[p](t)+i[pm](t))*(lambda[4](t)-lambda[3](t))+a*s(t)*(i[p](t)+i[pm](t))*(lambda[1](t)-lambda[2](t)))/(n.w[1]): u[2] := min(max(0, C), 1): C := (b*(i[m](t)+i[pm](t))*s(t)*(lambda[3](t)-lambda[1](t))+b*(i[m](t)+i[pm](t))*i[p](t)*(lambda[4](t)-lambda[2](t)))/(n.w[2]): u[3] := min(max(0, j), 1): j := (i[p](t)*lambda[2](t)+i[pm](t)*(lambda[4](t)-lambda[7](t))-(i[p](t)+i[pm](t))*lambda[5](t))/w[3]: u[4] := min(max(0, o), 1): o := (i[m](t)*lambda[3](t)-i[pm](t)*(lambda[4](t)-lambda[7](t))-(i[m](t)+i[pm](t))*lambda[6](t))/w[4]: u[2] := 0: u[3] := 0:

sys := diff(s(t), t) = Q+delta[1]*r[p](t)+delta[2]*r[m](t)+delta[3]*r[pm](t)-(a*(1-u[1])*(i[p](t)+i[pm](t))/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+mu)*s(t), diff(i[p](t), t) = (1-u[1])*a*(i[p](t)+i[pm](t))*s(t)/n-(1-u[2])*b*(i[m](t)+i[pm](t))*i[p](t)/n-(sigma[1]+u[3])*i[p](t)-(alpha[1]+mu)*i[p](t), diff(i[m](t), t) = (1-u[2])*b*(i[m](t)+i[pm](t))*s(t)/n-(1-u[1])*a*(i[p](t)+i[pm](t))*i[m](t)/n-(sigma[2]+u[4])*i[m](t)-(alpha[2]+mu)*i[m](t), diff(i[pm](t), t) = (1-u[2])*b*(i[m](t)+i[pm](t))*i[p](t)/n+(1-u[1])*a*(i[p](t)+i[pm](t))*i[m](t)/n-(y+u[3]+u[4])*i[pm](t)-(alpha[1]+alpha[2]+mu)*i[pm](t), diff(r[p](t), t) = (sigma[1]+u[3])*i[p](t)+(e*y+u[3])*i[pm](t)-(delta[1]+mu)*r[p](t), diff(r[m](t), t) = (sigma[2]+u[4])*i[m](t)+(y*g*(1-e)+u[4])*i[pm](t)-(delta[2]+mu)*r[m](t), diff(r[pm](t), t) = (y*(1-g)*(1-e)+u[3]+u[4])*i[pm](t)-(delta[3]+mu)*r[pm](t), diff(lambda[1](t), t) = lambda[1](t)*(a*(1-u[1])*(i[p](t)+i[pm](t))/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+mu)-lambda[2](t)*(1-u[1])*a*(i[p](t)+i[pm](t))/n-lambda[3](t)*(1-u[2])*b*(i[m](t)+i[pm](t))/n, diff(lambda[2](t), t) = -c[1]+lambda[1](t)*(1-u[1])*a*s(t)/n-lambda[2](t)*((1-u[1])*a*s(t)/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+sigma[1]+u[3]+alpha[1]+mu)-lambda[4](t)*(b*(1-u[2])*(i[m](t)+i[pm](t))/n+(1-u[1])*a*i[m](t)/n)-lambda[5](t)*(sigma[1]+u[3]), diff(lambda[3](t), t) = -c[2]+lambda[1](t)*(1-u[2])*b*s(t)/n+lambda[2](t)*(1-u[2])*b*i[p](t)/n-lambda[3](t)*(u[4]+alpha[2]+sigma[2]+a*(1-u[1])*(i[p](t)+i[pm](t))/n-(1-u[2])*b*s(t)/n)-lambda[4](t)*((1-u[2])*b*i[p](t)/n+a*(1-u[1])*(i[p](t)+i[pm](t))/n)-lambda[6](t)*(sigma[2]+u[4]), diff(lambda[4](t), t) = -c[3]+lambda[1](t)*((1-u[1])*a*s(t)/n+(1-u[2])*b*s(t)/n)-lambda[2](t)*((1-u[2])*b*i[p](t)/n+(1-u[1])*a*s(t)/n)-lambda[2](t)*((1-u[2])*a*i[m](t)/n+b*(1-u[1])*s(t)/n)-lambda[4](t)*(y+u[3]+u[4]+alpha[1]+alpha[2]+mu-(1-u[2])*b*i[p](t)/n-b*(1-u[1])*i[m](t)/n)-lambda[5](t)*(e*y+u[3])-lambda[6](t)*(y*g*(1-e)+u[4])-lambda[7](t)*(y*(1-g)*(1-e)+u[3]+u[4]), diff(lambda[5](t), t) = -lambda[1](t)*delta[1]+lambda[5](t)*(delta[1]+mu), diff(lambda[6](t), t) = -lambda[1](t)*delta[2]+lambda[6](t)*(delta[2]+mu), diff(lambda[7](t), t) = -lambda[1](t)*delta[3]+lambda[7](t)*(delta[3]+mu), i[p](0) = 300, i[m](0) = 200, i[pm](0) = 150, r[p](0) = 200, r[m](0) = 150, r[pm](0) = 150, s(0) = 1000, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0, lambda[4](T) = 0, lambda[5](T) = 0, lambda[6](T) = 0, lambda[7](T) = 0:
p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = 0.1e-3,maxmesh=2400);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging
 

Rewrite your equation as  k*x = 3/2*arctan(x)  and look at 2 graphs  y = k*x  and  y = 3/2*arctan(x) . Since the derivative of the second function at x=0 is 3/2, we get the obvious answer:   k∈(-∞,0]∪[3/2,+∞) . Below you can see the plot for k=3/2 :

 plot([3/2*x, 3/2*arctan(x)], x=-5..5, -5..5, color=[red,blue]);  

Edit. The solution above is true for the case b = 0. For arbitrary b<>0, we rewrite the equation  f(x)=b  as  k*x - b = 3/2*arctan(x) . Using the properties of functions and geomeric considerations, we get the final answer (k0(b) is the value of  k  for which for a given  b  the graphs touch each other ):

For  -3*Pi/4<b<0  or  0<b<3*Pi/4 :   k∈(-∞, 0] ∪ [k0(b), +∞) 
For  b<=-3*Pi/4  or  b>=3*Pi/4 :   k∈(-∞, 0) ∪ (0, +∞)

Below we can see the animation for the case b = -1. Note that the value of k0(b) can only be found numerically.

restart:
fsolve({k0*x+1=3*arctan(x)*(1/2),k0=3/2*1/(1+x^2)});
plots:-animate(plot,[[-3*Pi/4,3*Pi/4,3*arctan(x)*(1/2),[t*cos(s),1+t*sin(s),t=-5..5]], x=-5..5,-5..5,color=[black$2,blue,red], thickness=[1$2,2$2], linestyle=[3$2,1$2]], s=arctan(eval(k0,%))..Pi, frames=100);
                  {k0 = 0.3313385522, x = 1.878055290}