A:=d*n0/dt+d*epsilon*n1/dt+d*epsilon^(3/2)*n2/dt+d*epsilon^2*n3/dt+epsilon*d*n0*u1/dx+epsilon^(3/2)*d*n0*u2/dx+d*n0*epsilon^2*u3/dx+d*epsilon^2*n1*u1/dx+d*epsilon^(5/2)*n1*u2/dx:
coeff(A,epsilon^(3/2));

                                     d*n2/dt+d*n0*u2/dx        

The list  L2 contains all the options for representing numbers from your list in a sum of no more than 5 different squares.

For some reason, the document is not displayed inline. Therefore, I copied the code and the output fragment:

restart;
Squares:=[seq(i^2,i=0..12)]:
k:=0:
for n from 129 to 129+13^2 do
for m from 1 to 5 do
for c in combinat:-choose(Squares,m) do
if n=add(c) then k:=k+1; L[k]:=[n,c] fi;
od: od: od:
L:=convert(L,list):
L1:={seq(p[1]=`+`(seq(sqrt(p[2][i])^`2`,i=1..nops(p[2]))),p=L)}:
L2:=[ListTools:-Categorize((x,y)->lhs(x)=lhs(y),[L1[]])]:
for p in L2 do
lhs(p[1])=`or`(seq(rhs(p[k]),k=1..nops(p)));
od;

Download sums_squares2.mw

Edit.

I made 2 graphs: the first is static with a cutout in the first octant for greater clarity, the second is dynamic (animation). The upper and lower surfaces are made in different colors and the scale is equalized in all axes (scaling=constrained):

restart;
f:=x^(1/2):
g:=x^2/8:
X:=r*cos(phi): Y:=r*sin(phi):
plot3d(eval([[X,Y,f],[X,Y,g]],x=r), r=0..4, phi=Pi/2..2*Pi, color=["LightBlue","Khaki"], scaling=constrained, axes=normal, labels=[z,x,y], orientation=[25,75], lightmodel=light1);
plots:-animate(plot3d,[eval([[X,Y,f],[X,Y,g]],x=r), r=0..4, phi=Pi/2..a, color=["LightBlue","Khaki"], scaling=constrained, axes=normal, labels=[z,x,y], orientation=[25,75], lightmodel=light1], a=Pi/2..5*Pi/2, frames=60);

You mean the arithmetic mean? Why waste your time and look for the appropriate command if the problem is solved simply by adding the distances between all kinds of pairs of vertices and dividing by the number of these pairs. Here is an example:

with(GraphTheory):
G := Graph([1, 2, 3, 4, 5], {{1, 2}, {1, 3}, {1, 4}, {1, 5}, {2, 3}, {2, 5}, {3, 4}, {4, 5}}):
Dists:=AllPairsDistance(G);
`+`(seq(seq(Dists[i,j],j=i+1..5),i=1..4))/10;

In this example (when we set the graph  G ), distances are not specified and, by default, the length of each edge is assumed to be 1.

Edit. The last line of code can be written shorter:
  add(add(Dists[i,j], j=i+1..5), i=1..4)/10;

Index:=proc(v::{Vector,Array})
local i0,i1,i;
i0:=`if`(v::Vector,1,lhs(op(v)[1]));
i1:=`if`(v::Vector,op(v)[1],rhs(op(v)[1]));
for i from i0 to i1 do
if v[i]>=0 then return i fi;
od;
end proc: 

Examples of use:
V:=Vector([-1,-1,-2,0,3,4]):
Index(V);
V1:=Array(0..3, [-1,-2,4,-5]):
Index(V1);
V2:=Vector([-1,-1,-2]):
Index(V2);
                                 4
                                 2

Edit.
                                       

f:=x->-(x^3-6*x^2+11*x-6)*(x^2-4*x+5)/(x^4-8*x^3+24*x^2-32*x+17);
plot([f(x),2-x], x=-2..6, linestyle=[1,3],color=[red,blue], thickness=[2,0], scaling=constrained);

Addition. There is a rational function with polynomials of lower degrees satisfying all conditions:

f:=x->(-x^3+6*x^2-11*x+6)/(x^2-4*x+5);
convert(f(x),parfrac);
plot([f(x),2-x], x=-4..8, -5..5,  linestyle=[1,3],color=[red,blue], thickness=[2,0], scaling=constrained);

Here is a procedure for this:

Grouping:=proc(L::listlist)
local L1;
uses ListTools;
L1:=[Categorize((x,y)->x[2]=y[2], L)];
map(p->[p[1,2],map(l->l[1],p)], L1);
end proc:

Example of use:

MyNestedList := [   [g1, [1,2,3]]     ,   [g2, [0,1,3]]   ,  [g3, [1,2,3]]  ,  [g4, [0,1,3]]]:
Grouping(MyNestedList);
                   [[[1, 2, 3], [g1, g3]], [[0, 1, 3], [g2, g4]]]

I replaced  D  with  Daw  because  D  is the differential operator in Maple. We have a transcendental equation  m = 0.9058763835e-3*sqrt(6)/(x^2*m^2*exp(-26.00000000*x*m))-1/2-0.2884615385e-1/(x*m)  that we can solve only numerically. On the plot we see that for each  x   there are 2 values of  m . In the code we find  m  with the smallest modulus:
 

Download Equation.mw

restart;
F:=()->RandomTools:-Generate(choose([1$2, 2$5, 3, 4, 5])):
# Example of use
seq(F(), i=1..100);
ListTools:-Collect([%]);

4, 1, 4, 2, 1, 2, 1, 1, 4, 2, 2, 2, 5, 1, 2, 5, 2, 3, 2, 2, 2, 2, 

  1, 1, 2, 4, 2, 5, 1, 1, 2, 3, 5, 2, 1, 2, 2, 4, 1, 1, 1, 1, 4, 

  2, 2, 5, 2, 2, 1, 3, 2, 1, 2, 2, 2, 5, 2, 4, 2, 5, 2, 1, 2, 2, 

  4, 1, 5, 4, 2, 1, 3, 2, 3, 3, 5, 2, 2, 1, 2, 1, 2, 1, 2, 2, 3, 

  2, 2, 1, 5, 3, 5, 5, 5, 2, 2, 1, 2, 4, 2, 1
          [[1, 25], [2, 44], [3, 8], [4, 10], [5, 13]]
 

See this a toy example:

A:=a+b*x=0:
coeff(A, x);
coeff(lhs(A), x);

                       Error, unable to compute coeff
                                               b

 

Use  fsolve  with  initialpoint  option :

fsolve([Asubs/Bsubs = 1, 2*Esubs = h*`J__ω`[49]], {a = 2*10^(-20) .. 4*10^(-20), k = 4*10^8 .. 8*10^8});

                            {a = 3.31129407052*10^(-20), k = 6.47179093038*10^8}

I got this this initialpoint  {a = 2*10^(-20) .. 4*10^(-20), k = 4*10^8 .. 8*10^8}  from the implicitplot. If you want to automate this process for a large number of applications, then try using the  DirectSearch package, which you can freely download from Maple Application Center.

Here is a procedure that does this not only for lists, but also for expressions similar to polynomials. Formal parameters: P - a list or an expression of type `+`,  var is a list of symbols for which the procedure isolates nonlinear members.

IsolateNonlinearTerms:=proc(P::{list,`+`,polynom},var::list)
uses ListTools;
selectremove(p->degree(p)>1 or degree(p)<0, `if`(P::list,FlattenOnce(P),[op(P)]),var);
end proc:

Examples of use:

w:=[[z, y, x, x^3, 1, -5*y], [x*z, x*y, y, 2*x, 1], [x*z, 1/z^2, z, x*y]];
IsolateNonlinearTerms(w, [x,y,z]);
IsolateNonlinearTerms(x^3-3*x*y^2+y^4+x-5*y+4*z-10, [x,y]);
                         

restart;
pde1 := (y+z)*(diff(u(x, y, z), x))+(z+x)*(diff(u(x, y, z), y))+(x+y)*(diff(u(x, y, z), z)) = 0;
pdsolve(pde1, u(x,y,z), HINT=strip);

See help on  pdsolve  command  for details.

Edit. Below are 2 more options with animation. In the first case, we plot the function of 2 variables for each value of c1. The third coordinate of the highest point gives you the desired maximum. In the second option, we directly plot the maximum as a function of c1:

Download question_on_animate_with_background_new2.mw

Download question_on_animate_with_background_new1.mw

Download plot5.mw