Let me explain more about what I’m trying to do.You can see the book “optical properties of photonic crystals” written by K.Sakoda (or any other books related to photonic crystals and plane wave expansion method) on this subject. I have an eigenvalue equation sum(kappa(G-Gprime)*(abs(K+Gprime)^2)*E(Gprime),Gprime)-omega^2*E(G)=0 (2.61) G and Gprime are 2D vectors of the form (m,n) m,n can be any integer. K is also a 2D vector but in a positive range K=(k1,k2). Omega is the frequency . I want to find omega versus k1 and k2 for circular dielectric rods placed in a square lattice. For such a structure kappa is of the form kappa(G<>0)=2*f*(1/epsilon[a] -1/epsilon[b])*BesselJ(norm(G)*alpha,1)/norm(G)*alpha (2.98) kappa(G=0)=f/epsilon[a] +(f-1)/epsilon[b] (2.99) norm(G)=sqrt((g1)^2+(g2)^2) G=(g1,g2) f, epsilon[a],[b], alpha are some numeric constants. The coefficients of E(G) in (2.61) make a matrix. The number of Gprime vectors make the dimension of the matrix. After obtaining the matrix as I was able to obtain, I should put its determinant to zero ,solve it for omega and then plot omega versus k1 or k2 to plot a curve called dispersion curves in physics like Fig 2.4 in the book I mentioed. By increasing the number of vectos(the dimension of the matrix) I will have better and more accurate curves. I hope these details would be enough. Equations and the diagram fig.2.4

I raised the precision (evalf[50](…),evalf[100](…)) but unfortunately the results for sol1 and sol2 still have a large difference (especially for larger amounts of a, b, c) . About your suggestions on plotting: I can only put one of the variables(e.g.k2) to numeric values, then I should solve the equation det(M)=0 for omega. After the olve command I will obtain omega1(k1), omega2(k1),..,omegan(k1) (n=the number of solutions). I think what you've suggested is something like a pointplot; instantiating both k1 and k2 before evaluating det will give me points which are not related to each other and I can't categorize them as omega1(k1), omega2(k1),.. . I don't know how to obtain reasonable curves out of such a pointplot. A part of Maple codes and the point plot obtained for a 7*7 matrix.

I'm brave enough to examine every method that might help me. Your idea works well for a 5*5 matrix but for larger matrices it doesn't seem well enough. For a 7*7 matrix, radnormal doesn't give a zero result (or even near zero). Idea 2 results in evalf(…sol1,..); -0.0002583…. evalf(…sol2,..); -0.0000119…. Idea 3 for Determinant(evalf[50]( eval(m,[a=0.05,b=0.03,c=0.07]) )); yields -0.0000120 Besides I don't understand what the code Normalizer:=x->x really means. To me it acts like a miracle, it reduces the computation time for LUDecomposition very much but unfortunately as I mentioned above, it doesn't yield the desired result. About the original problem I can say I want to evaluate the determinant, put it to zero for I could be able to solve it for c and then plot real solutions for c vercus a or b or both( for more information on the last part you can see my comment “please help me..” on the post “how to extract real solutions” .)

It's a dense matrix in 3 variables a,b,c and the entries of the matrix can be for example (0.1702*Pi-1)*sqrt(abs(a-1)^2+abs(b-1)^2)-13*c^2 or -0.12*Pi*sqrt(a^2+b^2+2*b+1) The dimension of the matrix is 13*13 or more. The coefficient of c^2 is the rank of the matrix (here 13) and just it appears in the diagonal part of the matrix. Also I have the assumptions assume(a>=0); assume(b>=0); Does such a matrix need such a large memory you've indicated in the post "memory allocation" ? Thanks Sara

Would you please explain more about LU decomposition?

I will try the codes you've written and study on it. Thanks a lot for your help

I will try the codes you've written and study on it. Thanks a lot for your help

What do you mean by quantifier elimination?

What do you mean by quantifier elimination?

Thanks a lot for your advice; it helped me a lot. Using 'piecewise' or 'if clause' doesn't matter, both gives the same result. But using 'vectors' instead of 'list' is the leading point. I had tried vectors for this problem before but I had trouble with 'piece..' or 'if cla..' ;they didn't understand a condition like " if x='a vector' " ,they just understand "if x= 'a number' ". So after your comment I worked more on it and finally I could solve my problem as follows: >with(LinearAlgebra): >len:=v→VectorNorm(v,2): >kappa:=y→piecewise(len(y)=0, 5, 1/len(y)): >Q[0]:= <0|0>: Q[1]:= <0|1>: Q[2]:=<2|0>: >S:=proc(R::Vector) local n,T,U,W,j; global kappa; for n from 0 to 2 do T[n]:=VectorAdd(R, −Q[n]); U[n]:=kappa(T[n]); W[n]:=U[n]*e(Q[n]); od; sum(W[j], j=0..2); end; >S(Q[1]); 2e([0 0])+5e([0 1])+2√5/5e([2 0]) No error message. Thanks for your help. (Q[n] are written as row vectors but here they aren't typed well,)

I've tried add. it doesn't work here.