We can easily find all matrices that commute with A by simply reducing the problem to solving a system of 4 equations with 4 unknowns:

restart;
A := Matrix(2, 2, [[3, 1], [-3, 2]]):
C := Matrix(2, 2, [[x, y], [u, v]]);
A.C=~C.A;
Sys:=convert(%,set);
solve(Sys);
C:=eval(C, %); # The general formula for C,where  x and y are any numbers
E:=eval(C, [x=1,y=2]); # The specific example
E.A=A.E;  # Check

The new procedure  NewLegend1  allows you to place the opaque rectangle anywhere on the plot. Please note that in the list  G , you must specify a spacing for the variable  y  also as for  x .

restart;
NewLegend1:=proc(f::algebraic,g::algebraic,R::list,G::list)
local A1, A2, A3, A4, B, a, b, l1, l2, T;
uses plottools, plots;
A1:=implicitplot([y=f,y=g], subsop(1=(x=op([2,1],G[1])..R[1,1]),2=(y=op([2,1],G[2])..R[1,2]),G)[], labels=["x","y"]);
A2:=implicitplot([y=f,y=g], subsop(1=(x=R[1,1]..op([2,2],G[1])),2=(y=op([2,1],G[2])..R[2,2]),G)[], labels=["x","y"]);
A3:=implicitplot([y=f,y=g], subsop(1=(x=R[2,1]..op([2,2],G[1])),2=(y=R[2,2]..op([2,2],G[2])),G)[], labels=["x","y"]);
A4:=implicitplot([y=f,y=g], subsop(1=(x=op([2,1],G[1])..R[2,1]),2=(y=R[1,2]..op([2,2],G[2])),G)[], labels=["x","y"]);
B:=rectangle(R[1],R[2], style=line);
a:=R[2,1]-R[1,1]; b:=R[1,2]-R[2,2];
l1:=line(R[1]+[0.1*a,-0.3*b],R[1]+[0.4*a,-0.3*b],color=rhs(G[3])[1],thickness=rhs(G[4]));
l2:=line(R[1]+[0.1*a,-0.7*b],R[1]+[0.4*a,-0.7*b],color=rhs(G[3])[2],thickness=rhs(G[4]));
T:=textplot([[(R[1]+[0.75*a,-0.3*b])[],y=f],[(R[1]+[0.75*a,-0.7*b])[],y=g]], font=[times,bold,14]);
display(A1,A2,A3,A4,`if`(nops(R)=2,B,NULL),l1,l2,T);
end proc:

Example of use:

NewLegend1(sin(x),cos(x), [[0.2,0.9],[1.4,0.3]], [x=0..2*Pi, y=-1..1, color=[red,blue], thickness=2, size=[900,400], scaling=constrained, axes=box]);

The  NumberTheory:-Divisors  command returns the set of all divisors of an integer. Procedure below solves your problem:

restart;
a:=proc(n)
local S, N, i, j, k;
uses NumberTheory;
k:=0;
S:=NumberTheory:-Divisors(n); N:=nops(S);
for i from 1 to N-2 do
if ormap(j->irem(S[j],S[i])=0, [$i+1..N-1]) then k:=k+1; fi;
od:
k;
end proc:

seq(a(n), n=1..24);

                               0, 0, 0, 1, 0, 1, 0, 2, 1, 1, 0, 3, 0, 1, 1, 3, 0, 3, 0, 3, 1, 1, 0, 5

Edit. The procedure  a  can be written much shorter:

a:=n->nops(NumberTheory:-Divisors(n))-nops(NumberTheory:-PrimeFactors(n)) - 1:

The third way: the procedure  a  can be written without calling the  NumberTheory  package, only using the kernel command  ifactors :

a:=proc(n)
local L, N;
L:=ifactors(n)[2];
N:=nops(L);
mul(L[i,2]+1, i=1..N)-N-1;
end proc:

Below is a possible variant of the desired procedure. The list  R  defines the position of the rectangle with legends. If the third element of the list  R  is some symbol, for example 'off', then the rectangle is not drawn.

restart;
NewLegend:=proc(f::algebraic,g::algebraic,R::list,G::list)
local A, B, a, b, l1, l2, T;
uses plottools, plots;
A:=plot([f,g], G[], labels=["x","y"]);
B:=rectangle(R[1],R[2], style=line);
a:=R[2,1]-R[1,1]; b:=R[1,2]-R[2,2];
l1:=line(R[1]+[0.1*a,-0.3*b],R[1]+[0.4*a,-0.3*b],color=rhs(G[2])[1],thickness=rhs(G[3]));
l2:=line(R[1]+[0.1*a,-0.7*b],R[1]+[0.4*a,-0.7*b],color=rhs(G[2])[2],thickness=rhs(G[3]));
T:=textplot([[(R[1]+[0.75*a,-0.3*b])[],y=f],[(R[1]+[0.75*a,-0.7*b])[],y=g]]);
display(A,`if`(nops(R)=2,B,NULL),l1,l2,T);
end proc:

Examples of use:

NewLegend(sin(x),cos(x),[[3.5,0.9],[4.3,0.3]],[x=0..2*Pi,color=[red,blue],thickness=2,size=[900,400],scaling=constrained, axes=box]);

NewLegend(exp(x),x,[[2.5,1.8],[3.6,1],'off'],[x=-2..5,color=[red,blue],thickness=2,size=[500,500],scaling=constrained, view=[-2..5,-2..5],axes=normal]);

Edit. 1. I don't understand why you color the rectangle and place it in the top left corner every time. If you want the procedure to be universal, you cannot specify in advance the place where the rectangle should be located. It may turn out that just in this place there will be something interesting on the plot, and the rectangle will block it. Also, with your way, you won't be able to use the  view  option which is often useful (see your second example with  exp(x) ).
2. The procedure is written for the case of two graphs, but it is easy to modify it to suit any number of graphs, for example, if you have only one function or 3 functions, etc.
3. Please note that all the options required for plotting are included in the list  G . I understand your desire to reduce the number of formal parameters, but the degree of automation should be reasonable. The main thing in my opinion is the best end result.

To extract single entries or entire rows or columns use square brackets instead of parantheses:

restart;
A:=Matrix([[1,2],[3,4]]);
A[1,..];
whattype(%);
A[..,1];
whattype(%);

For rows  A[1,..] = A[1]  and so on.

It looks like

restart;
f:=(x,t)->sqrt(2)*exp(-0.6*(0.5*x+0.01*t)^2):
plot([f(x,0),f(x,20),f(x,40)], x=-10..10, color=[yellow,black,blue], thickness=3, size=[900,500]);

What 2 lines did you mean? In the plot from your worksheet, I see only one straight line and 2 coordinate axes.
General advice: to find the intersection point of two lines  F(x,y) = 0  and  G(x,y) = 0 , solve the system

solve({ F(x,y) = 0,  G(x,y) = 0}, {x,y});

Edit. Sorry, I didn't notice your second plot right away. So do as Carl wrote or just

restart;
solve({y=5*x-1, y=-x+7});

restart; 
P := 2^n-1: 
select(isprime, [seq(eval(P, n = k), k = 1 .. 100)]);

  [3, 7, 31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111]

Edit. If you do not need prime numbers of the form  P  themselves, but exponents  n , then below is another code:

restart; 
P := 2^n-1: 
select(p->isprime(eval(P, n=p)), [$1..100]);

                 [2, 3, 5, 7, 13, 17, 19, 31, 61, 89]

restart;
a:=3: R:=7: z:=sqrt(R^2-(r-a+R)^2):
plot3d([[r,phi,z],[r,phi,-z]], phi=0..2*Pi, r=2..3, coords=cylindrical, scaling=constrained);

The formula  sqrt(R^2-(r-a+R)^2)  was obtained by elimination  phi  from the system  { z=R*cos(phi), r=a-R*(1-sin(phi)) } 

Edit. I added the option  scaling=constrained so that you can watch the form change depending on the parameters change.

In my opinion the simplest way is to use a triple for-loop:

restart;
C := y^2*z + y*z^2 = x^2:
for x in {0,1} do
for y in {0,1} do
for z in {0,1} do
if C then print(['x','y','z']=[x,y,z]) fi;
od: od: od:

To plot a 3-variable equation, you must use  plots:-implicitplot3d  command. To get the best quality we can first express  z  through  x  and  y and then plot:

restart;
eq:= x^2-y^2*z-y*z^2;
plots:-implicitplot3d(eq,x=0..1,y=0..1,z=-7..7, style=surface, grid=[40,40,40]);
solve(eq, z);
plot3d([%], x=0..1,y=0..1, style=surface);

Edit. A more interesting picture is obtained if we increase the ranges for  x  and  y :  x=-1..1, y=-1..1. It is worth noting one drawback of  plot3d  command in this example: Maple does not draw a surface near z-axis. This is due to the fact that for small values of  y , the value of  z  tends to infinity. The use of  plots:-implicitplot3d  command does not have this drawback. We get good quality when used together  plot3d  and  plots:-implicitplot3d :

restart;
eq:= x^2-y^2*z-y*z^2;
A:=plots:-implicitplot3d(eq,x=-1..1,y=-1..0.01,z=-5..5, style=surface, grid=[50,100,50]):
Sol:=solve(eq, z);
B:=plot3d([Sol], x=-1..1,y=-1..1, view=-5..5, style=surface, axes=normal, grid=[500,500]):
plots:-display(A,B, view=[-1..1,-1..1.2,-5..5]);

Let's help Maple a little. If  t^2 - 4 = s^2 - 4  and  t<>s  then  s = -t . Below we get that there will be infinitely many double points:

restart;
x:=t->sin(t): y:=t->t^2-4:
solve({sin(t)=sin(s),s=-t}, AllSolutions);
eval(%, _Z1=k);
[x,y]=eval([x(t),y(t)], %) assuming k::posint;  # All the double points
# k = 1, 2, 3, ...

                               [x, y] = [0, Pi^2*k^2-4]

The numerical solution:

restart;
x:=t->sin(t): y:=t->t^2-4:
Sys:=[x(t)=x(s),y(t)=y(s)]:
plots:-implicitplot(Sys, t=-10..10,s=-10..10, color=[red,blue], gridrefine=3);
fsolve(Sys, {t=-4..-1,s=1..4});
fnormal(eval([x(t),y(t)], %));

                                    [0., 5.869604404]

Because there are many solutions, then the plot is needed to isolate the root.

With local, you made a new name (it's not clear why?) from the command  pdetest . Just replace the line local pdetest := (ans, de1)  with the line  test := pdetest(ans, de1)

restart;
did:=D^14+23*D^13+144*D^12-30*D^11;
subsindets(did,symbol^integer,t->(D@@degree(t))(y)(x));
convert(%, diff);

You can use a common  if ... fi:  for this:

restart;
str:="A";
x:=10;
if x=10 then str:=cat(str," it was 10"); x:=11 else
str:=cat(str," it was not 10"); x:=9 fi;

Edit.