In fact, if  sin(t)=a  and  cos(t)=b, then from this does not follow that  t=arctan(a/b)  (only tan(t)=a/b). The correct equality will be  t=arctan(a,b) (of course up to a summand a multiple of 2 * Pi). See help on two-argument arctan  ?arctan . 
For proof, it is sufficient to use the command  expand  and no assumptions are required:

Download Q20201108_new.mw

Edit: another way:

Matrix(6,2, (i,j)->`if`(j=1,rand(10..20)(),rand(50..100)()));

Download RM.mw

You can use  the empty symbol  ``  for this  (to prevent parentheses from expanding):

restart;
eq:=a*``(2*A + 2*B) + b*C = c*``(2*A + B) + d*``(B*2) + e*``(C + 4*B);
subs({a = 6, b = 1, c = 6, d = 1, e = 1}, eq);

Your system has an infinite number of solutions depending on one parameter  C2 :

Download sys_1.mw

The procedure  PS  finds the first  m  members of this sequence.

PS:=proc(m)
local L, S, k, n, a, b;
L:=[]; 
for k from 1 to m do
for n from 1 do
a:=ithprime(n); b:=cat(op(L),a);
if not (a in L) and isprime(parse(b)) then
L:=[op(L),a];  break fi;
od; od;
op(L);
end proc: 

Example of use:
PS(100);
2, 3, 11, 7, 41, 31, 17, 163, 23, 79, 197, 241, 29, 37, 59, 193, 227, 229, 239, 439, 929, 337, 257, 1447, 509, 19, 293, 1723, 1619, 937, 179, 367, 251, 1063, 4241, 1291, 521, 1951, 443, 139, 191, 1753, 1217, 673, 53, 883, 809, 109, 5381, 3733, 311, 967, 449, 2683, 9239, 577, 971, 1453, 101, 1051, 107, 3301, 233, 1627, 2729, 2689, 5669, 1597, 2153, 1117, 5099, 10429, 5303, 6991, 4253, 1777, 281, 349, 5273, 1543, 1997, 1867, 1109, 769, 167, 4261, 3821, 14923, 173, 1531, 71, 13267, 14699, 3877, 137, 2389, 3803, 2671, 5879, 1039

Edit.
 

restart;
pde := -e*(diff(u(x, y), x, x) + diff(u(x, y), y, y) )+ a*diff(u(x, y), x) - x*(-p^2*x + 1)*sin(l*y) = 0;
A := -p; B := (1-2*a*A)/p; C := (2*e*A - a*B); p := e*l^2;
eval(pde, u(x,y) = (A*x^2 + B*x +C)*sin(l*y));
simplify(%);

Download check.mw

In addition to the numerical optimization commands  Optimization:-Maximize  and  Optimization:-Minimize, Maple has symbolic optimization commands  maximize  and  minimize , but none of these commands work with functions containing parameters. The simple examples below illustrate this. Optimization:-Maximize  and  Optimization:-Minimize commands does not work with piecewise-functions, even if they do not contain parameters.

Download optimization_examples.mw

Check if maximize  and  minimize commands will work with your function for specific parameter values. If so, you can write a procedure to calculate these, the parameters of this procedure will be the parameters of your function. This procedure will allow you to solve a variety of problems, for example, to plot graphs the dependence of the maximum and minimum on parameters.

Let  a  and  b  be the major and minor semiaxes of the ellipse, c - focal distance. From the condition  b=с  and the equality  a^2=b^2+c^2 , we easily get that  a=c*sqrt(2)  and the eccentricity  e=1/sqrt(2) .  By the condition b = c  the ellipse is determined up to similarity, the alpha angle does not change, therefore, in the calculations, we can assume that  c=1 .
Below  r1  and  r2  are the distances from a point  M(x,y)  on the ellipse to the ellipse foci.

Download max_angle.mw

Use  op(eqs)  or  eqs[ ]  for this.

restart;
S:={a*x,b*x, c*(x+y)};
S1:=op~(S);
S1 minus {a,b,c};

restart;
eq:=sqrt(y)=tanh(x);
eq^2;

Edit.

Download PMSM_eq_new1.mw

In each specific example, the shortest way is probably to use the  subsop  command:

restart;	
sol:={v[1]=t,v[2]=3/2*t,v[3]=v[3]}:
subsop(3=(),sol);

R:=1.33*10^(-9)*C/m^2;
evalf[3](op(1,R)*``(`*`(op(2..3,R))));

restart;
expr:=exp(x)^3*sin(x)+3/(exp(x)^n):
A:=exp(t::anything)^n::anything=exp(t*n):
applyrule([1/A,A], expr);

Edit.