@mehdibaghaee 
It is your turn to say something about the points 1-2 above. After that I will be able to comment. (How do you know that the method works without examples with known eigenvalues?)

@mehdibaghaee 

The problem is changed now.
You should test the problem with small matrices. Nobody will want to see (and work with) your 700x700 matrices.
Note that Eigenvalues (unlike Determinant) works for big matrices (due to special algorithms) so probably the change of variables is not needed.

Download lista1.mw

@sand15 

assuming positive  means that all the variables are assumed to be >0.
On my system (Maple 2016.2, 64 bit) I get FAIL  for     is( r^2 - r*s + s^2> 0) assuming r >0, s>0;

positive means real and >0. So,  is(I, positive)  is of course false (not being real).

r^2+s^2 > 0   is simple for Maple (sum of two positive numbers); s^2 - s*r + r^2  seems to be too complicated for the moment.

Edit (concerning "simple conditions"):

is( (r - s)^2 + r^2  + s^2  > 0)  assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive, r>s/2;
                              true

@mehdibaghaee 

My question was about the type of the problem. Where/how  those 10^6 variables are they used?
For instance, do you have a system of equations with 10^6 unknowns?

@mehdibaghaee 

It would be interesting to know the type of the problem which needs 10^6 variables and can be solved.

@Markiyan Hirnyk 

But you have spent a rather long time with this problem as I can see. So, selectively busy.

@Markiyan Hirnyk
It is obvious (as I have already said) that the procedure cannot print an infinity of points.
And I do not have enough spare time to analyse more deeply the problem (for which there are many articles in the literature).
But you are free (and invited) to write a better procedure. I will vote it up!

@Zeineb 

Using distributions (such as Dirac) is a delicate task. They are not actually functions but Maple allows algebraic operations with them. It is the user's responsibility to ckeck that these operations have a mathematical sense (for example, usually it is not possible to multiply a non-regular distribution by a function which is not C^oo).

For t=7Pi/6 AND 11Pi/6

@Markiyan Hirnyk 

@Markiyan Hirnyk 

Why don't you try to understand the problem first?
t runs in the interval [0,2Pi).
For example, the point [2,0] is double (in this approach): it is obtained for t=0 AND for t=Pi.

@Markiyan Hirnyk 

In a Lissajous curve m,n should be relatively prime i.e. igcd(m,n)=1. Otherwise all the points are double/multiple and the procedure shows only the "extreme" ones. So, use instead dblpoints(1,2,2,3).

For a non-trivial procedure see the chapter "Advanced examples" in

Char B.W. et al - First Leaves, A Tutorial Introduction to Maple V, Springer 1992

It is very fast and works for huge numbers.

@Carl Love 

Hello Carl, glad to "see" you here again.
I think that OP's confusion is due to the fact that t was omitted in 'maximize'.
So, probably it is better to write:

f:= cos(2*t/m) + cos(2*(t+5)/m):
plot('maximize'(f,t), m= 1..10);

@Christopher2222 

It is not the loop. If you recompute several type evalf(add(...)) ==> random results.
And it is not the add, the same thing happens for `+`(seq(...)).
The problem is the randomness, not necessarily the 3 digits!
This is not nice.