It sounds like the PDEs equivalent of 'DEtools[odeadvisor]'. If so, I am not aware of something like that.

As GAMMA(a)-GAMMA(a,1)=Int(exp(-t)*t^(a-1),t=0..1) (check ?GAMMA),
abs(GAMMA(a)-GAMMA(a,1))<=Int(t^(a-1),t=0..1)=1/a->0 as a->infinity.

With Maple, you can do some steps as:

J:=convert(GAMMA(a)-GAMMA(a,1),Int) assuming a>0;
IntegrationTools:-Combine(subs(infinity=inf,J));
                      1
                     /
                    |                (-1 + a)
                    |   exp(-_k1) _k1         d_k1
                    |
                   /
                     0

('Combine' seems to have problems with 'infinity')

is(abs(exp(-_k1))<=1) assuming  _k1>=0,_k1<=1;
                                 true

value(Int(_k1^(-1+a),_k1 = 0 .. 1)) assuming a>0;
	1/a

was not aware of a thing named "Legendre economization". Does it mean that the truncated Taylor series, expressed in terms of these polynomials, is computationally more economical?

may be what you are looking for.

PS. However, it sounds like the output of

surfdata(data2, axes = boxed);

is not the intended surface. This seems better:

PLOT3D(POLYGONS(data2),STYLE(PATCHNOGRID),AXESSTYLE(BOX),
ORIENTATION(6.,68.));

I do not find a statement in ?pdsolve,numeric that the number of independent variables is limited by 2 (if missing, it should be added), but there are many hints about it.

It produces an infinite recursion (phi calls phi...):

phi :=   (x1,x4) -> phi(x1,x4):
phi(x1,x4);

Error, (in phi) too many levels of recursion

I prefer GRTensor for this kind of calculations.

Here, Maple assumes implicitly that R, omega, etc, are complex. So you have to assume them real and positive:

z := R*omega^2*C^2/(omega^2*C^2+R^2)+I*(R^2*omega*C
+omega^3*C^2*L+R*omega^2*C^2)/(omega^2*C^2+R^2);

abs(z) assuming positive;
  /   2      4  4        2                3  2            2  2 2\1/2
  |  R  omega  C       (R  omega C + omega  C  L + R omega  C ) |
  |----------------- + -----------------------------------------|
  |      2  2    2 2                     2  2    2 2            |
  \(omega  C  + R )                (omega  C  + R )             /


argument(z) assuming positive;
                    2                3  2            2  2
                   R  omega C + omega  C  L + R omega  C
            arctan(--------------------------------------)
                                       2  2
                                R omega  C

this 'command':

with(PDEtools):
ODE:=12*diff(f1(x),x$5)^(7/2)*diff(f2(x),x$4)^9
+12*diff(f3(x),x$3)*diff(f1(x),x$2)^9:

command:=proc(ode,f)
indets(ode,specop(anything,diff));
select(has,%,f);
max(op(map(difforder,%)));
end proc:

command(ODE,f1);
                                  5

command(ODE,f2);

                                  4

command(ODE,f3);

                                  3

It does not work because it is not a subexpression of 'ex1', but collect for them does work. So, "constructively" you can do:

ex1 := C10*C20*a^2+2*C10*C20*a*b+C10*C20*b^2-m*u^2*a*C10
+m*u^2*b*C20+1/2*e1*e2*a^2+e1*e2*a*b+1/2*e1*e2*b^2:
collect(ex1,[C10,C20,e1,e2],distributed,factor):
subs(a+b=z,%):
collect(%,z,factor):
subs(z=a+b,%);

          /          e2 e1\        2      2
          |C20 C10 + -----| (a + b)  - m u  (a C10 - C20 b)
          \            2  /

You can do this:

expand(f,3*t);

                  sin(3 t) cos(5) - cos(3 t) sin(5)

You may trace what 'SumTools:-DefiniteSum:-Frontend' does this way:

kernelopts(opaquemodules=false):
infolevel[all]:=5:
trace(SumTools:-DefiniteSum:-Frontend):

you can do:

sort(data,(a,b)->a[3] < b[3]);

                  [[3, 6, -1], [2, 3, 1], [5, 0, 3]]

solves dgl.

It can be done this way:

F:= -1/2*eta^3+3/2*eta;
plots:-contourplot(x*F,x=-2..2,eta=-1..1);

May be you mean cylindrical coordinate lines, something that can be approached this way:

with(plots):with(plottools):
P:=changecoords(plot3d(r^2,r=0..1,phi=0..2*Pi,axes=normal),cylindrical):
C:=coordplot(polar,labelling =true,scaling=constrained):
f := transform((x, y) -> [x, y, 0]):
display({P, f(C)}, orientation=[46, 52]);