> ListTools:-Flatten([[1,15,21],[5,7,33]]);

And then, as Georgios suggests, you can use sort on the result.

Yes.  Working in floating point to a certain number of digits, Maple can only produce a result with that many digits of precision.  It's doing the best that can be done under those circumstances.  In your example, with the default Digits = 10:

> fsolve(tan(x) = x, x = 4 .. 6);

      4.493409458

The actual solution is somewhere between 4.493409457 and 4.493409458, but no numbers strictly between 4.493409457 and 4.493409458 can be represented with Digits = 10.

tan(4.493409458) is 4.493409460, a bit too high, while tan(4.493409457) is 4.493409439, a bit too low.
You might notice that the fact that the derivative of tan(x) - x at the root is rather high (over 20) makes it rather unlikely that you would get a solution where Maple's values for tan(x) and x are exactly the same: tan(x) - x will change by about 20 ulps from one float representable with a certain number of digits to the next one.

Yes.  Working in floating point to a certain number of digits, Maple can only produce a result with that many digits of precision.  It's doing the best that can be done under those circumstances.  In your example, with the default Digits = 10:

> fsolve(tan(x) = x, x = 4 .. 6);

      4.493409458

The actual solution is somewhere between 4.493409457 and 4.493409458, but no numbers strictly between 4.493409457 and 4.493409458 can be represented with Digits = 10.

tan(4.493409458) is 4.493409460, a bit too high, while tan(4.493409457) is 4.493409439, a bit too low.
You might notice that the fact that the derivative of tan(x) - x at the root is rather high (over 20) makes it rather unlikely that you would get a solution where Maple's values for tan(x) and x are exactly the same: tan(x) - x will change by about 20 ulps from one float representable with a certain number of digits to the next one.

It looks like a bug in match.  In perhaps its simplest form:

> match(f/(f+1)^2 = f/(f-d)^2, f, 'S');

                    false

I'm submitting an SCR.

It looks like a bug in match.  In perhaps its simplest form:

> match(f/(f+1)^2 = f/(f-d)^2, f, 'S');

                    false

I'm submitting an SCR.

See the replies to CKOne in www.mapleprimes.com/forum/testingprimes

See the replies to CKOne in www.mapleprimes.com/forum/testingprimes

It's rather hard to tell what the shape of the surface should be from such few points.  But the same fitting procedures (either LSSolve or Statistics[Fit]) will work for any type of surface, specified as z = f(x,y)  for some function f with certain parameters (preferably occurring linearly).  So for example in Erik's solution you might replace p1 * x1 + p2 * x2 + p3 by p1 * x1^2 + p2 * x1 * x2 + p3 * x2^2 + p4 * x1 + p5 * x2 + p6.

> H := Statistics[Fit](p1 * x1^2 + p2 * x1 * x2 + p3 * x2^2 + p4 * x1 
   + p5 * x2 + p6, independentDataMatrix, LIPO, [x1, x2]);

It's rather hard to tell what the shape of the surface should be from such few points.  But the same fitting procedures (either LSSolve or Statistics[Fit]) will work for any type of surface, specified as z = f(x,y)  for some function f with certain parameters (preferably occurring linearly).  So for example in Erik's solution you might replace p1 * x1 + p2 * x2 + p3 by p1 * x1^2 + p2 * x1 * x2 + p3 * x2^2 + p4 * x1 + p5 * x2 + p6.

> H := Statistics[Fit](p1 * x1^2 + p2 * x1 * x2 + p3 * x2^2 + p4 * x1 
   + p5 * x2 + p6, independentDataMatrix, LIPO, [x1, x2]);

Another way, if you are just interested in the minimum and maximum, is to use member to find the position of the minimum or maximum value in a list of values.  This is then the position of an element in the original list that produces the minimum or maximum value.  Thus:

> X := [3,4,-1,2,4]:
  f := x -> x^2:
  fX  := map(f, X);  
  maxf:= max(fX);
  member(maxf, fX, 'jmax'):
  xmax:= X[jmax];

Another way, if you are just interested in the minimum and maximum, is to use member to find the position of the minimum or maximum value in a list of values.  This is then the position of an element in the original list that produces the minimum or maximum value.  Thus:

> X := [3,4,-1,2,4]:
  f := x -> x^2:
  fX  := map(f, X);  
  maxf:= max(fX);
  member(maxf, fX, 'jmax'):
  xmax:= X[jmax];

Like this:

> m := 20;
  n := 20;
  Total := 0;
  for o from 0 to 20 do
    a:= (factorial(m)*factorial(n)*factorial(o)*(sum(sum(sum(sum(sum(sum(f, p =
0 .. limp), q = 0 .. limq), r = 0 .. limr), s = 0 .. lims), t = 0 ..
limt), u = 0 .. limu))^2/2^(m+n+o));
    Total := Total + a;
   end do;

Like this:

> m := 20;
  n := 20;
  Total := 0;
  for o from 0 to 20 do
    a:= (factorial(m)*factorial(n)*factorial(o)*(sum(sum(sum(sum(sum(sum(f, p =
0 .. limp), q = 0 .. limq), r = 0 .. limr), s = 0 .. lims), t = 0 ..
limt), u = 0 .. limu))^2/2^(m+n+o));
    Total := Total + a;
   end do;

> Abcds:= [seq([comb[4*n-3],comb[4*n-2],comb[4*n-1],comb[4*n]], 
     n=1..N/4)];
  Error := map(Err1 @ op, Abcds);
  Emin := min[defined](Error);
  member(Emin,Error,'p');
  abcdmin:= Abcds[p];

            abcdmin := [.1, 0, 1, 0]

> Abcds:= [seq([comb[4*n-3],comb[4*n-2],comb[4*n-1],comb[4*n]], 
     n=1..N/4)];
  Error := map(Err1 @ op, Abcds);
  Emin := min[defined](Error);
  member(Emin,Error,'p');
  abcdmin:= Abcds[p];

            abcdmin := [.1, 0, 1, 0]