Question 1:
I'm not sure there is a limitation of the number of plots you can display in a single plot3d command. For instance

N := 100:  # try 1000 if you want
F := [seq(i, i=1..N)]:
r := rand(0. .. 1.):
c := () -> ColorTools:-Color([r(), r(), r()]):
plot3d(F, x=0..1, y=0..1, color=[seq(c(), n=1..N)])

I would say the limitation is more related to human capabilities than to software's. Are you capable to visually separate let's say, 10 surfaces drawn together, or to undersand what they represent ?
Personally I begin having problems with 3 surfaces.

Question 2:
I suggest you to run one of these commands and see what kind of graphics Maple can do

plots[interactive]();

# or 
plots[interactive](abs(solnum));

There is an obvious difference at line 5: in aa1.mw Num:=100 whilie Bum:=200 in aa2.mw

Other differences:
Open aa1.mw, from the menubar: File > Export As > Maple Text (produces a file bamed aa1.txt).
Do the same for aa2.mw

Use some online file or text comparator if you don't have one, for instance Compare (in this case open the txt file and drag-and-drop the texts in the appropriate regions)
The differences are blue highlighted (aa1 is on the left)



So, contrary to what you write, these two programs are NOT identical

Maybe "identical programs" mean for you "the two programs, while different, should give the same result" ? 
But this is another question. 

So check he parts of the codes where differences occur and make yourself sure that those parts give the same results.
For instance are you sure these definitions of T11 have similar effects, see aa.mw

Here is a revised version of your aa2 code where I switch from aa1 code to aa2 code for T11, ..., T16 
This is why you get different results : aa2_mmcdara.mw
It is obvious that T11 (at least) has different expressions in the two codes (look after the plot):

# "T11 from aa1" = "T11 from aa2" iif k0 = -I, which is never the case:

k0 := (.5*Pi+(2*(N-1))*Pi/(Num-1))/(l1+2*l2+2*l3+2*l4)

            0.29089 + 0.00186 (N - 1) Pi

So, contrary to what you write, these two processes are NOT similar

Solve {E2, E3, E4} wrt {b, c, d} (for instance).

Plug the solution into the equality E1 - f(Ea) = 0.

Build an objective function J(T)(A) depending on the remaining parameter a (noted A), equal to the integral of  (E1 - f(Ea))²  over t from 0 to some value T (or over any other range).

Minimize J(T)(A) wrt to A for any given value of T (J(T)(A) represents the L2 norm of E1 - f(Ea) for the value a=A).

Here is an example

Download Aph1_mmcdara.mw

Download Q_KKT_mmcdara_1.mw

restart
N := 3:
S := [T, P]:
mul(mul(cat~(S, `__`, ||i)), i=1..N);

# or
# mul(mul(cat~(S, `_`, ||i)), i=1..N)

result:


mul.mw

If you really want the symbol "_*" to appear:
 

Download mulP.mw

Two stories:

1. In French universities, Maple homeworks usually go something like this: the teacher gives the students a problem, and they have two options: either stay at the university after class and use Maple on the university's workstations, or go home and write the code in  text format (or in Maple if they've bought a student license). They then deliver their work in mw or txt format.
   After what "someone" corrects the job and tell them the possible error they produced.
   And if the student gets errors he has to correct them before delivering another code.
   A friend of mine who teaches Statistics at Bordeaux University told me that quite often students go on specific forums to ask for users to help them fix their code instead of trying to do that b themselves.
2. It happens this "someone" it is a thesis student which does the correcting job to earn a few money. This same friend of mine told me about the case of one of his student who asked on a forum to fix issues on R homework codes that it was intended to correct himself but was too lazzy to do so.

Are you one ofthese two kinds of guys : the student waiting for a "une bonne poire" (I allow myself this french term you surely understand... for others I found "sucker" or "soft touch" using slang translators) to correct his paper, or the lazzy teacher assistant ?
I'm pretty sure the fact you stubbornly refuse for years to upload any Maple worksheet has something to do with that.

Finally, concerning your question, is this what you want?

 

If it is so, here are some hints for correction

• First of all start by cleaning up your code (twice the same in a row).
   
• Make a choice: either you want to use the geometry package, either you don't.
   
• You wrie a procedure named distance but distance is also a function of the geometry do if you chosed to use geometry do not write a procedure named distance.
   
• You are mixing notations several times: for instance circle(I1, r1) is the  plottools definition of a circle, not a geometry definition.
   
• geometry has functions incircle to find the incircle of a given triangle and excircle to find the find three excircles of a given triangle: why don't simply use them (all the more that your procedure plot_exscribed_circles builds wrong excircles).
   
• For ratio (R) = 1 your "Apollonius hyperbola" is a straight line whatever A and B:

  f := (x, y) -> sqrt((x - A[1])^2 + (y - A[2])^2)/sqrt((x - B[1])^2 + (y - B[2])^2) - 1:
  collect~(isolate(f(x, y)=0, y), x)
  
                                      2       2       2       2
           (2 A[1] - 2 B[1]) x   -A[1]  - A[2]  + B[1]  + B[2] 
     y = - ------------------- - ------------------------------
             2 (A[2] - B[2])            2 (A[2] - B[2])        
  

  For any other positive value of R it is a circle:

  f0 := sqrt((x - A[1])^2 + (y - A[2])^2)/sqrt((x - B[1])^2 + (y - B[2])^2) = R;
                                  (1/2)          
                         / 2    2\               
                         \x  + y /               
                   -------------------------- = R
                                        (1/2)    
                   /       2          2\         
                   \(x - 4)  + (y - 2) /         
  
  f0 *~ denom(lhs(f0)): 
  %^~2: 
  (lhs-rhs)(%): 
  tx := add(select(has, [op(expand( %))], x)):
  ty := add(select(has, [op(expand(%%))], y)): 
  sx := Student:-Precalculus:-CompleteSquare(tx, x):
  sy := Student:-Precalculus:-CompleteSquare(ty, y): 
  
  add(map(s -> s/(-R^2+1), [op(sx+sy)])): 
  f1 := map(add, (select=-remove)(has, [op(%)], {x, y}))
  
                         2                2             
            /        2  \    /        2  \          4   
            |     4 R   |    |     2 R   |      20 R    
            |x + -------|  + |y + -------|  = ----------
            |      2    |    |      2    |             2
            \    -R  + 1/    \    -R  + 1/    /  2    \ 
                                              \-R  + 1/ 
  

Example of what you could/should do:
(and I consider I'm already "une bonne poire" doing this) here is way (not the simplest one!) to draw the incircle of triangle ABC

plot_inscribed_circle := proc(__A, __B, __C)
  local A, B, C, T, Ii:
  point(A, op(__A)):
  point(B, op(__B)):
  point(C, op(__C)):
  triangle(T, [A, B, C]);
  incircle(Ii, T,'centername'=o);
  draw(Ii, style = line, color = green, thickness = 2)
end proc:

plot_inscribed_circle(A, B, C)

Now it's up to you to do some effort... or not.

Stop sending us a mess of c... like you just did and telling us "hey guys, deal with this".
Have a little respect for people whose passion is helping others.

If any moderator thinks I've crossed a red line by writing this reply, may it feel free to delete it. 

Do not write 

# Equation of the line passing through A and B
    lineAB := sort(Equation(line(k, [A, B], [x, y])));

but

# Equation of the line passing through A and B
    line(k, [A, B], [x, y]):
    lineAB := Equation(k);

instead (do not embed the definition of the line within Equation)

Tou will stiill get an error for item=[1, 6, 12] because points A and B are confounded (so line k is not defined).

item;
detail(A);
                           [1, 6, 12]
           GeometryDetail(["name of the object", A], 

             ["form of the object", point2d], 

             ["coordinates of the point", [-2, -8]])

 detail(B);
           GeometryDetail(["name of the object", B], 

             ["form of the object", point2d], 

             ["coordinates of the point", [-2, -8]])

It is up to you to fix this problem

Here is an idea to avoid the error when line k is undefined 
 

Download Corrections_mmcdara.mw

Nevertheless this is a step-by-step Maple code which mimics what you would do by hand (note the result is exact and no complex expression ever occur... which is a necessary condition for the result to be found by hand). My first version: Another_way.mw

The version above uses solve, this new one is even simpler as it doesn't:
 

Download Another_way_variant.mw

A first code:  P is made of 15 lists, S is the sum of the elements of each list:

L := [a, b, c, d]:
P := combinat:-choose(L):
P := subsop(1=NULL, P):
S := map(add, P):
print~(S):

                               a
                               b
                             a + b
                               c
                             a + c
                             b + c
                           a + b + c
                               d
                             a + d
                             b + d
                           a + b + d
                             c + d
                           a + c + d
                           b + c + d
                         a + b + c + d

Permuting all the lists in P gives a new collection PP of lists. But adding the elements of permuted lists always give the same result.
For instance add([a, b]) = a+b  and  add([b, a]) = a+b.
If you want the addition of the elements in [b, a] gives b+a you must cheat, for instance:

PP := map(op@combinat:-permute, P):

f := proc(p):
  if numelems(p) = 1 then
    cat(`#mo("`, p[1], `")`)
  else
    cat(`#mrow(`, seq( cat(`mo("`, op(i, p), `"),mo("+"),`), i=1..numelems(p)-1 ), `mo("`, op(-1, p), `"))`)
  end if:
end proc:

# example (partial result)

PP[1..4];
print~(f~(PP[1..4])):
                   [[a], [b], [a, b], [b, a]]

All results here  permutations.mw
 

@JAMET 

There are 3 types of conics depending on the sign of the determant of matrix Q:

Q := < <A | B>, <B | C>>

# for a quadratic form written  A*x^2 + 2*B*x*y + C*y^2 + ...

In each type there exist several "degenerated" situations (for instance a parabola can reduce to two parallel straight lines, or an ellipse to  a single point).
Generally the classic rotation+translation transformation applies if ysed with precaution..

These special cases (cf Bibmath et Wiki for more details about special cases) make delicate the writting of a code aimed to build the reduced form in all situations.

For your conic:

(
This code Gender_and_Signature.mw determines the type and the signature of your quadratic form:

• Determinant       ==> type = parabola,
• Signature=(2, 1)  ==> sub type = parabola non degenerated parabola)

)

Reduced form and siplays for your conic
 

Download Reduced_Parabola_2.mw



Ice on the cake

This code returns the type and subtype of a quadratic form (Note: I'm not sure about these designations; maybe one should use "class" instead and dropt types and subtypes ?)
 

Download Conic_Classification.mw

UPDATED 09/14/2024

This is typically a problem which is not well posed as it enables an infinity of solutions.

To make you understand this let's take this mechanical analogy

You are in 1D. A rod is made of two rods of same length L, the lewt one weights 10 kg, the right one 1kg.
The left extremity of the left rod is at position x=1 (so the system extends from x=1 to x=2).
Compute the position of the center of mass (the mean in Statistics) ?
Would you say: that it is (10*L/2 +1*(L+L/2)) / 11 = 13*L/2
if so you would be wrong: you missed a fundamental assumtion about the density of the mass distribution in each rod. 
The answer 13*L/2 is correct iif you assume the density is constant in each bar.

Relation with your problem.
If the probability density function (coninuous case) or the probability mass function (discrete case) is constant over each interval, then you can resume each interval by its "center of mass" anf rewrite your problem in termes of "weighted data": 

 

So your problem has an infinity of answers unless you correctly formulate assumptions on what happens within each class, or on the distribution of the underlying population those classes are built from.


More importantly: I hope the question you pose is a simple amusement because in real data analysis you never compute statistics from data organized in classes.
You compute them from the raw data that you can organized further on as classes... if you want.
Computing statistics on classes and numbers implies loss of information and ambiguity.
This loss of information prevents from calculating statistics for class-based data.


Nevertheless here is my contribution to your problem stat2.mw.

I UPDATED this reply by adding file BOUNDS.mw which provides bounds for the mean, variance and quartiles when no assumption are made on the distribution of the data or when you do not know them:that's the only objective answer an honorable person should give you.
Whatever tthose assumptions might be your are assured that those bounds are these ones:

_____________________________________________________________________________________

Auxilliary point concerning your MMA results:

At the evidence the person who answered your question reduced you problem to a problem of assessing statistics on weighted data (which is wrong).

For your information, R is one of the World reference code in statistics; its Weighted.Desc.Stat package is aimed to compute statistics for weighted data. 
Considering data consist in the vector x centers of classes  (which corresponds to the homogenous mass density distribution in the mechanical analogy) and in vector mu of weights, on gets these estimators when x is the vector of the centers of classes:.

> x <- c(1.25, 3.75, 6.25, 8.75)
> mu <- c(18, 11, 13, 6)
> w.mean(x, mu)
[1] 4.114583
> w.var(x, mu)
[1] 7.028537
> w.sd(x, mu)
[1] 2.651139
> 

Note thatthe variance here differs from MMA's (=7.549370660)

The formula Weighted.Desc.Stat uses to assess this variance, and why it uses it, are reproduced below.

(1)  Replace all occurrence of 

myfac[i]

by

myfac(i)

(2)  Replace all occurrence sum by add.
      (look to the terms addition and sum in the help pages, Page types=Dictionary)

sum_vs_add.mw

Orignal reply deleted by the author.

Revised version: Updated.mw

 

From the reference you give:

 

Download Dihedral_Angle.mw

As you've already noticed, there's nothing immediate..
The trick to derive equations 9, 11 and 13 is to multiply both sides of equation 7 by G' and integrate wrt eta (without without forgetting the "constant" K)

(Done with Maple 2015)
 

Download This_way.mw