I understood the  "identify"  word as the need to collect identical objects into separate groups (lists). The code below does this:

restart;
Eq1:=5.09295817894067*`τu`[1, 1]-30.5577490736439*`τu`[2, 1]+178.253536262923*`τu`[3, 1]-30.5577490736439*`τu`[1, 2]+183.346494441862*`τu`[2, 2]-1069.52121757753*`τu`[3, 2]+178.253536262923*`τu`[1, 3]-1069.52121757753*`τu`[2, 3]+6238.87376920228*`τu`[3, 3]:
Eq2:=5.09295817894067*`τu`[1, 1]+10.1859163578814*`τu`[2, 1]+15.2788745368241*`τu`[3, 1]-30.5577490736439*`τu`[1, 2]-61.1154981472883*`τu`[2, 2]-91.6732472209439*`τu`[3, 2]+178.253536262923*`τu`[1, 3]+356.507072525849*`τu`[2, 3]+534.760608788841*`τu`[3, 3]-3/7:
Eq3:=5.09295817894067*`τu`[1, 1]-30.5577490736439*`τu`[2, 1]+178.253536262923*`τu`[3, 1]+10.1859163578814*`τu`[1, 2]-61.1154981472883*`τu`[2, 2]+356.507072525849*`τu`[3, 2]+15.2788745368241*`τu`[1, 3]-91.6732472209439*`τu`[2, 3]+534.760608788841*`τu`[3, 3]-9/7:
Eq4:=5.09295817894067*`τu`[1, 1]+10.1859163578814*`τu`[2, 1]+15.2788745368241*`τu`[3, 1]+10.1859163578814*`τu`[1, 2]+20.3718327157631*`τu`[2, 2]+30.5577490736484*`τu`[3, 2]+15.2788745368241*`τu`[1, 3]+30.5577490736484*`τu`[2, 3]+45.8366236104784*`τu`[3, 3]-12/7:
Eq5:=5.09295817894067*`τu`[1, 1]-30.5577490736439*`τu`[2, 1]+178.253536262923*`τu`[3, 1]+50.9295817894067*`τu`[1, 2]-305.577490736439*`τu`[2, 2]+1782.53536262923*`τu`[3, 2]+504.202859715131*`τu`[1, 3]-3025.21715829077*`τu`[2, 3]+17647.1000900295*`τu`[3, 3]-18/7:
Eq6:=5.09295817894067*`τu`[1, 1]+10.1859163578814*`τu`[2, 1]+15.2788745368241*`τu`[3, 1]+50.9295817894067*`τu`[1, 2]+101.859163578814*`τu`[2, 2]+152.788745368241*`τu`[3, 2]+504.202859715131*`τu`[1, 3]+1008.40571943027*`τu`[2, 3]+1512.60857914560*`τu`[3, 3]-3:
Eq7:=5.09295817894067*`τu`[1, 1]+10.1859163578814*`τu`[2, 1]+15.2788745368241*`τu`[3, 1]-30.5577490736439*`τu`[1, 2]-61.1154981472883*`τu`[2, 2]-91.6732472209439*`τu`[3, 2]+178.253536262923*`τu`[1, 3]+356.507072525849*`τu`[2, 3]+534.760608788841*`τu`[3, 3]-3/7:
Eq8:=5.09295817894067*`τu`[1, 1]+10.1859163578814*`τu`[2, 1]+15.2788745368241*`τu`[3, 1]+10.1859163578814*`τu`[1, 2]+20.3718327157631*`τu`[2, 2]+30.5577490736484*`τu`[3, 2]+15.2788745368241*`τu`[1, 3]+30.5577490736484*`τu`[2, 3]+45.8366236104784*`τu`[3, 3]-12/7:
Eq9:=41.7622570673196*`τu`[3, 1]+41.7622570673196*`τu`[1, 3]+15.2788745368220*`τu`[1, 1]+83.5245141346398*`τu`[2, 3]+30.5577490736443*`τu`[2, 1]+113.063671572516*`τu`[3, 3]+83.5245141346398*`τu`[3, 2]+30.5577490736443*`τu`[1, 2]+61.1154981472892*`τu`[2, 2]:

S:=[seq([i,Eq||i], i=1..9)];
S1:=[ListTools:-Categorize((u,v)->u[2]=v[2], S)];
map(t->op~(1,t), S1);
op(select(t->nops(t)>1, %));

The final result:     
                                    [2, 7], [4, 8]

If you just need to remove duplicates, then do

{seq(Eq||i, i=1..9)};

If we divide the total area of all the necessary pieces of plywood into the area of 1 standard plywood sheet  48*96=4608, we get about 2.56. Thus, at least 3 standard sheets are required. I just cut out proportional sized pieces of paper and in 10 minutes I found the right solution. Thus, 3 standard sheets are needed. All calculations and visualization are below. For ease of identification, pieces of the same size are painted in the same color. All sizes are in inches.

Plywood1.mw  

Use the  isolve  command for this:

Addition. Unfortunately,  isolve  command is not effective enough and quite often (in more difficult situations) it does not cope with the task of finding integer points. On this occasion, look at these 2 posts:  https://www.mapleprimes.com/posts/202437-Integer-Points-In-Polyhedral-Regions  and  https://www.mapleprimes.com/posts/202542-Integer-Points-In-Curved-Regions

Download help_list_integer_equation_new.mw

This is the standard integer linear programming problem. See wiki  https://en.wikipedia.org/wiki/Cutting_stock_problem 

See also  How do I solve Cutting Stock Problem in Maple - MaplePrimes 

We project points  C  and  D  on the plane of the base of the cylinder. We get points  C'  and  D' . Obviously  ABC'D'  is a rectangle. For the origin we take the center of this rectangle, we direct the  x axis through the middle of  AB, and the  y axis through the middle of  BC' . Let  AB = a  and  C'B=b .  To find  a  and  b , we have a simple system of two equations with 2 unknowns. We solve it with Maple:

solve({a^2+b^2=6^2, b^2+2^2=a^2, a>0, b>0});

               {a = 2*5^(1/2),  b = 4}

Further finding the coordinates of all points is obvious.

If you want to enter those programmatically (without a palette directly from the keyboard) in 1D Input or 2D Input modes, then this can be done as in the examples below:

conjugate(a); 
`#mover(mi("a"),mo("~"))`;
`#mover(mi("a"),mo("ˆ"))`;
`#mover(mi("a"),mo("→"))`;

Output:
                                     

You can write even shorter with the same result:

restart;
f := (Pi*x+2*c+2*m)/(mu__c*S)+2*epsilon/(mu__a*S)+(Pi*x+2*c)/(mu__s*S);
int(1/f, x=0..w)   assuming positive;

You can simplify the final result a little:

combine(%) assuming positive;

I removed the double quotes and added 2 more options:

restart;
f:=x->1/x^2:
plot(f(x), x = -10 .. 10, thickness = 3, labels = ["", ""], tickmarks = [[4 = infinity, -4 = -infinity], [0]], view = [-6 .. 6, 0 .. 5], scaling = constrained, size = [900, 375]);

If you mean the  FunctionAdvisor  command then it appeared in Maple 9 in 2003. 
See  ?updates,v9

restart;
Sys:={x=a/(1+C^2)^(3/2), y=a*C^3/(1+C^2)^(3/2)}:
simplify(map(t->t^(2/3), Sys));
`+`(op(%));
applyop(simplify,2,%) assuming a>0, C>0;

In fact, of course, the original set  {x=..., y=...}  does not coincide with the set  x^(2/3)+y^(2/3)=a^(2/3) , since for a chosen parameter  a  sign of  x  will always coincide with the sign of  a  (in my interpretation without the minus). However, under the conditions  a>0, C>0, x>=0, y>=0  these sets coincide.

Try

f:=(l::list)-> eval([y,y*z-x,-15*x*y-x*z-x],zip(`=`,[x,y,z],l));

f([1,2,3]);
                              [2, 5, -34]

But why not just write

f:=(x,y,z)-> [y,y*z-x,-15*x*y-x*z-x];

f(1,2,3); 
                                 [ 2,5,-34]  

# Ellips x^2/a^2+y^2/b^2=1, a>b
# Circle x^2+y^2=a^2
# The ellipse is obtained by compressing the circle to the Ox axis with a coefficient k=b/a

a:=5: b:=3: k:=b/a:
F:=plottools:-transform((x,y)->[x,k*y]):
A:=plot([a*cos(t),a*sin(t), t=0..2*Pi], color=blue, thickness=2):
plots:-display(A, F(A));

Your problem is the problem of finding  the geometric median  for a finite set of points. See  https://en.wikipedia.org/wiki/Geometric_median

If not all points lie on one straight line, then the solution is always unique  as in your example . Otherwise (all points on one straight line) if there are an even number of points, then there are infinitely many solutions, as in the example below:

Points:=[0, 1, 3, 10]:
plot(add(abs(x-p), p=Points), x=-2..12);

Unfortunately, both commands below do not solve this simple example:

restart;
Points:=[0, 1, 3, 10]:
S:=add(abs(x-p), p=Points);
minimize(S, location);
Optimization:-Minimize(S);

Coordinates of points found wrong. You can take (for example) as below

with(geom3d):
point(B, -3/2, -2, 0):
point(A, 3/2, -2, 0):
point(C, -3/2, 2, 0):
point(S, x, y, z):
solve([distance(S, A) = 3, distance(S, B) = 5, distance(S, C) = 7], [x, y, z], explicit);

I just swapped the points  A  and  B .

The method you are referring to is based on the usual shift and has a significant drawback - the thickness of such a surface in different places will be different. This is clearly seen in 2D-examples below. In the first plot, the second curve is obtained by shifting the parabola  y=x^2  by 0.1 up. In the second graph, the second curve passes through the ends of the all normals with the length 0.1 at the points of the parabola  y=x^2 :

restart;
plot([x^2, x^2+0.1], x=-1.5..1.5, color=black, scaling=constrained);
A:=plot(x^2, x=-1.5..1.5, color=black, scaling=constrained):
B:=plot([x-2*x/sqrt(1+4*x^2)*0.1, x^2+1/sqrt(1+4*x^2)*0.1, x=-1.5..1.5], color=black, scaling=constrained):
plots:-display(A, B);

Below is the same method with normals applied to the surface defined by the explicit equation  z=x^2-y^2. For clarity, both sides of the surface and its end face are made in different colors:

f:=(x,y)->x^2-y^2:
A:=plot3d(f(x,y), x=-1..1, y=-1..1, color=pink, style=surface):
Dx:=diff(f(x,y),x): Dy:=diff(f(x,y),y):
d:=sqrt(Dx^2+Dy^2+1): h:=0.1:
B:=plot3d([x-Dx/d*h,y-Dy/d*h,f(x,y)+h/d], x=-1..1, y=-1..1, color=blue, style=surface):
C1:=plot3d([1-2/sqrt(2^2+Dy^2+1)*h*s,y-Dy/sqrt(2^2+Dy^2+1)*h*s,f(1,y)+h*s/sqrt(2^2+Dy^2+1)], s=0..1, y=-1..1, color=cyan, style=surface):
C2:=plot3d([-1+2/sqrt(2^2+Dy^2+1)*h*s,y-Dy/sqrt(2^2+Dy^2+1)*h*s,f(-1,y)+h*s/sqrt(2^2+Dy^2+1)], s=0..1, y=-1..1, color=cyan, style=surface):
C3:=plot3d([x-Dx/sqrt(Dx^2+2^2+1)*h*s,1+2/sqrt(Dx^2+2^2+1)*h*s,f(x,1)+h*s/sqrt(Dx^2+2^2+1)], s=0..1, x=-1..1, color=cyan, style=surface):
C4:=plot3d([x-Dx/sqrt(Dx^2+2^2+1)*h*s,-1-2/sqrt(Dx^2+2^2+1)*h*s,f(x,-1)+h*s/sqrt(Dx^2+2^2+1)], s=0..1, x=-1..1, color=cyan, style=surface):
plots:-display(A, B, C1, C2, C3, C4, axes=none);

Exactly the same method can be used for a surface defined parametrically, only the formulas for the normals will be different. For a surface specified implicitly, some technical difficulties arise in the implementation of this method, but I can roughly imagine how they can be solved. Soon I am going to write a post with a procedure to automate all the steps in this plotting.

thick.mw