The Answers of John and Kitonum have given a minor symptomatic solution to your problem but have failed to mention that your code can be vastly simplified.

If f and h are both sets, then all you need is

f intersect h;

The following works if f or h are sets or lists or a mixture:

select(`in`, f, h);

The print command is very, very rarely a good way to get your results. Whenever you feel tempted to use print, you should try to think of another way, or ask a Question here.

You want 20 arrays each with 20 elements? Use a 20x20 Matrix. Fill it with a single call to Sample:

X:= Matrix((20$2), datatype= float[8]):
N:= Statistics:-RandomVariable(Normal(0,1)):
Statistics:-Sample(N, X):

Now the jth column of X can be accessed as X[..,j] and is equivalent to a 20-element Vector.

(I know that the above commands can be shortened. I wrote it the way I did for pedagogy.)

Since q[2] is 0, this is trivial in that the sum reduces to one term. But I'll assume that you want to do this in more general cases. Also, I'll assume that x = x[1] and y = x[2]; otherwise, the answer is, of course, 0.

q:= <1,0>:
add(add(diff(a*x[1]*(1-x[1])-b*x[1]*x[2], x[j], x[k])*q[j]*q[k], j= 1..2), k= 1..2);

What ODE system? These are PDEs, not ODEs. You need to use pdsolve

One procedure will handle both of your questions, of course. Here it is:

Normalize:= proc(L::list)
local f:= proc() option remember; 0 end proc, k:= 0, x;
     forget(f);
     for x in L do
          if f(x)=0 then
               k:= k+1;
               f(x):= k
          end if
     end do;
     f~(L)
end proc:

Normalize([1,3,1,3,2,2,4,4]);
      [1, 2, 1, 2, 3, 3, 4, 4]

A question for the experts: Why is the forget necessary in this procedure? It gave erroneous results without it. I used Maple 2016. Surely this must be a major bug!

If you use the "raw" animated GIF export from Maple, then the display time for each frame will be 0.1 seconds and, it'll play as a continuous loop. You can download a GIF animation editor called Easy GIF Animator. With this, you can set the display time for frames individually, and you can turn off the continuous looping.

plot([
     piecewise(s < 2, s, 2+2*sin((s-2)/2)),   #x-coord
     piecewise(s < 2, 0, 2*(1-cos((s-2)/2))), #y-coord
     s= 0..2+Pi
]);

You asked:

What is the meaning of the suffix "+0.I"? Does it mean that there are further decimal digits which are not displayed?

I call these "spurious imaginary parts". They occur when an answer is real but its computation uses complex numbers. Consider the cubic p:= x^3 - 3*x - 1. A plot will quickly reveal that it has three real roots. But solve(p, x) will show that complex numbers are used to calculate those roots. The spurious imaginary parts can be safely ignored if you have reason to believe that the solutions are real. The +0.*I parts can be removed with simplify(..., zero).

How do the first and third eigenvalues, which are equal, result in different eigenvectors? As per my understanding, equal eigenvalues should have equal corresponding eigenvectors.

No, that's only true in a rare situation called a defective matrix. Consider: What are the eigenvectors of the 2x2 identity matrix? How about the 2x2 zero matrix? The usual situation is that there are as many linearly independent eigenvectors as the multiplicity of the corresponding eigenvalue.

The earliest date recognized by this system is January 1, 1901. So, to get to January 1, 1900, we need to subtract or add 365.

Finance:-AdvanceDate("January 1, 1901", 23021 - 365);

     "January 12, 1963"

Finance:-DayCount ("January 1, 1901", "January 10, 1963") + 365;

     23019

Apparently, your day count was off by two days. Remember that 1900 was not a leap year.

Download sorted_plot.mw

Here's an example:

restart:
ODEs:= {
     diff(x(t),t)=sigma*(y(t)-x(t)),
     diff(y(t),t)=x(t)*(rho-z(t))-y(t),
     diff(z(t),t)=x(t)*y(t)-beta*z(t)
}:
ICs:= {x(0)=0, y(0)=0, z(0)=1}:
params:= {sigma= 10, beta= 8/3, rho= 28}:
n:= 600:  step:= 0.1:
A:= dsolve(
     eval(ODEs union ICs, params),
     numeric,
     output= Array([seq(step*k, k= 0..n)])
):
DesiredOrder:= [t, x(t), y(t), z(t)]:
Column:= proc(e) local p; member(e,A[1,1],p); p end proc:
M:= A[2,1][.., [seq(Column(e), e= DesiredOrder)]]:
ExportMatrix("C:/MyMatrix.csv", M, target= csv):    

Note that A[2,1] is the actual output data, and A[1,1] is just a Vector that shows the order by which dsolve arranged the columns in A[2,1].

Is that less confusing than the other answers that you found?

Your Question is ambiguous as to whether you want abs(max(Fm(x))) or max(abs(Fm(x))). Please clarify. In the code below, I assume the latter.

Since you're using Excel, I'm going to guess that you're using Windows. If not, the directions below need to be modified a little.

Directions for installing a Maple package:

1. Go to the Maple Applications Center and download the DirectSearch package
2. Create a folder somewhere on your computer to hold Maple packages. In my case, that folder is named "C:/Maple_packages".
3. From the zip folder, select the subfolder of your desired language: English or Russian.
4. Copy the files DirectSeach.mla and DirectSearch.help to the new folder. Then copy the whole zip folder to the new folder, just in case you ever need it (unlikely).
5. If you have a Maple initilization file, open it with Notepad. If you don't have one, create a file named "C:/Program Files/Maple 2016/lib/maple.ini" (this assumes that you have Maple 2016) and open it with Notepad.
6. Put in that file the line libname:= "C:/Maple_packages", libname; using the folder name you created in step 2. Save the file. If you already have a Maple initialization file, then it doesn't matter where in the file you include the above line. Make sure that you spell the second libname correctly. If you don't, you won't be able to use Maple at all until you fix this.
7. If you have Maple open, then close and re-open it. If you use Maple strictly through Excel, then I guess there's no need for this step. But we need to insure that your Excel is starting a fresh Maple session. So close Excel and reopen it.

The Maple command to find the maximum is

DirectSearch:-GlobalOptima(abs(Fm(x)), [x = a..b], maximize)[1];

To find the minimum, change maximize to minimize.

For example,

eq4:= dsolve({eq2, incs}, numeric, method= rosenbrock, range= 0..10, maxfun= 0);
plots:-odeplot(eq4, [t, y(t)], 9..10, numpoints= 10000);

Please don't include irrelevant plaintext output in your Questions. It makes them difficult to copy-and-paste.

Use IntegrationTools:-Expand.

The command

plots:-complexplot3d(f(z), z= a+b*I..c+d*I);

is equivalent to

plot3d(
     abs(f(x+I*y)), x= a..c, y= b..d,
     color= (argument(f(x+I*y))+Pi)/2/Pi, labels= [Re('z'), Im('z'), ""]
);

Thus, to achieve what you want, you can use

plot3d(
     abs(f(x+I*y)), x= a..c, y= b..d,
     color= (-argument(f(x+I*y))+Pi)/2/Pi, labels= [Re('z'), Im('z'), ""]
);

I wish that this could be done by manipulating the color option in complexplot3d, but I can't make it work.

Edit:I updated the plot3d commands based on an analysis of the actual code of `plots/complexplot3d`.