T[0]  := f:
T[10] := sqrt(omega) * 12^(1/4);
T[20] := evalindets(T[0], specfunc(anything, exp), exp @ proc(P) global T; radnormal(op(1, P) / T[10]) * T[10] end):
T[30] := [sqrt(omega) = Omega / (12^(1/4)), (omega) = (Omega / (12^(1/4)))^2];
T[40] := factor(expand(subs(T[30], T[20])));
T[50] := [exp(Omega) = beta[1], exp(I * Omega) = beta[2]];
T[60] := factor(radnormal(subs(T[50], T[40]))):
op(map(rhs = lhs, T[50])), Omega = solve(T[30][1], Omega), T[60];
subs(%);

_interface(showassumed = 0):
assume(x > 0, t > 0, k > 0):
pde := diff(u(x, t), t) = k * diff(u(x, t), x$2):
ic := u(x, 0) = 0:
bc := u(0, t) = t:
sol:= pdsolve({pde, ic, bc}, u(x, t)):

`@`(sort, combine, expand, evalc, map)(Re, [solve](5*x^3 - 5*x + 1, x));

T:=table([(25, 1) = -39, (16, 151) = 32, (33, 1) = -54, (1, 1) = 29, (13, 1) = 32, (31, 101) = -7, (6, 51) = -10, (11, 101) = -1, (28, 151) = -39, (18, 51) = -65, (4, 151) = 29, (8, 151) = -10, (23, 101) = 23, (34, 51) = -54, (40, 151) = 87, (36, 151) = -54, (9, 1) = -1, (37, 1) = 87, (21, 1) = 23, (14, 51) = 32, (22, 51) = 23, (20, 151) = -65, (27, 101) = -39, (3, 101) = 29, (19, 101) = -65, (24, 151) = 23, (32, 151) = -7, (30, 51) = -7, (38, 51) = 87, (7, 101) = -10, (10, 51) = -1, (29, 1) = -7, (35, 101) = -54, (17, 1) = -65, (26, 51) = -39, (15, 101) = 32, (12, 151) = -1, (39, 101) = 87, (5, 1) = -10, (2, 51) = 29]);

F := proc (T) Matrix((op @ map2)(`@`(min .. max, op, map2), op, [1, 2], ({indices})(T)), proc (i, j) if assigned(op(procname)[i, j]) then op(procname)[i, j] else 0 end if end proc[T]) end proc;

G := proc (T) Matrix((op @ map2)(`@`(min .. max, op, map2), op, [1, 2], ({indices})(T)), convert(op(op(T)),set)) end proc;

F(T);
G(T);

member((a-2*b-2*c)^2*(2*a-b+2*c)^2, M, 'i'), i;

A := log[12](27) = a, log[36](24) = b;
B := seq(simplify(e), e = A);
C := frontend(eliminate, [{B}, {b, ln(2)}], [{Non(function)}, {}]);
b = convert(subs(C[1], b), parfrac, a);

F := proc(S, T)
   local A, n, m, x, y, v, w, V, W, i, f, g;
   n := `$`(1 .. nops(S));
   m := `$`(1 .. nops(T));
   f := () -> args[1];
   g := `@`(assign, op, map);

   g(V @ f = op, [n], 'S');
   g(W @ f = op, [m], 'T');
   g(v = `@`(simplify, eval, V), [n]);
   g(w = `@`(simplify, eval, W), [m]);

   A := frontend([solve], [(x = 'v')~({n}) union (y = w)~({m}), indets(w~({m}), function) union y~({m})], [{Non(function)}, {}]);
   A := (x = V) ~ ([n]), subs(A[1], y~([m]));

   (W ~ ([m])) =~ (subs(A))
end;

F('{log[12](27)}, {log[36](24)}');
F('{log[2](3), log[3](5), log[7](2)}, {log[140](63)}');

plots[display](
plot3d(arctan(y, x), y = -1 .. 1, x = -1 .. 1, style = point, grid = [121, 121], symbol = POINT, color = black),
plots[display](plottools[cuboid](evalf([-1, -1, -Pi]), evalf([1, 1, Pi]), color = red, thickness = 3), style = wireframe),
plots[spacecurve]([0, x, arctan(0, x)], x = -1 .. 1, color = blue, thickness = 3, numpoints = 61), axes = boxed, scaling = unconstrained, orientation = [125, 35]);

judge:=proc(x,y,t) if not x = y then t := evalb(x < y) end if end proc;
W := proc(L1, L2)
   local i, t, X;
   if nops(L1) = nops(L2) then
      X := 'judge(op(i, L1), op(i, L2), 't')';
      for i while not assigned(t) do eval(X) end do;
      t
   else evalb(nops(L1) < nops(L2)) end if
end proc;
ops := m -> sort(op(2, OperandsTable(m)), (a, b) -> W([lhs](a), [lhs](b)));
ops(m);;

function_coeffs := proc(A, n::set(name))
local S, T, v;
   S := indets(A, {function});
   if nargs < 2 then v := {} else v := n; S := select(has, S, v); end if;
   T := {Non(map(identical, S))};
   frontend(proc(A, S) local V; [coeffs](collect(A, S, distributed), S, 'V'), [V] end proc, [A, S union v], [T, {}])
end proc;

eq2 := -4*A[2]*cos(2*x)-16*A[4]*cos(4*x)-36*A[6]*cos(6*x)-64*A[8]*cos(8*x)+a*A[0]+cos(8*x)*a*A[8]+cos(6*x)*a*A[6]+cos(4*x)*a*A[4]+2*cos(2*x)*q*A[0]+cos(2*x)*a*A[2]+q*A[8]*cos(6*x)+q*A[8]*cos(10*x)+q*A[6]*cos(4*x)+q*A[6]*cos(8*x)+q*A[4]*cos(2*x)+q*A[4]*cos(6*x)+q*A[2]*cos(4*x)+q*A[2];

zip(proc(C,F) if not hastype(F, trig) then C * F fi end, function_coeffs(eq2));

function_coeffs(sin(x)+cos(y)+x^2+x*y+y^2+2*x+1, {y});

typematch(Heaviside(x),'Heaviside(a::algebraic)');

epsilon:=0.01;
display(
seq(spacecurve([cot(phi)/sin(phi), 1.61, phi], phi = r, numpoints = 3000, axes = normal, coords = spherical, color = red, axes = boxed),r=(-Pi+epsilon .. -epsilon,epsilon..Pi-epsilon))
);

simplify([f0,f1,f2],{a*b},[a,b]);
#simplify(simplify([f0,f1,f2],{a^2-q1,b^2-q2,a*b},[a,b,q1,q2]),{a^2-q1,b^2-q2,a*b},[q2,q1,b,a]);
#subs(q1=a^2,q2=b^2,map(limit,simplify([f0,f1,f2],{a^2-q1,b^2-q2,a*b-r},[a,b,q1,q2,r]),r=0));

x||(1..4)||y||(1..4);
``||(a,b,c,d)||(1..4);

#old syntax
#x.(1..4).y.(1..4);
#``.(a,b,c,d).(1..4);