@vv Thanks for the trick! Much appreciated

@tomleslie Thanks.

The purpose of the code is to come up with an expression for Tcr as a function of (Aiso, Do, Gc, a, b, c, mu, t).
Note (Aiso, Do, Gc, a, b, c, mu, t) are known constants and Tcr is quad root (fourth root) and I need all the roots. For example, if Aiso=0 from that expression, you can obtain a desire solution in less than a second. But unfortunately Aiso is not zero

@JohnS You're absolutely right. The bc(s) are for the "simple support" conditions for an orthotropic plate bending problem.

I realized I can solve this problem using a variational principle approach without using a double Fourier expansion of w(x,y) for a 1D case (i.e. beam). And I'm trying to extend the same principle for the 2D case. But the problem now is how to apply my bs(c) in Maple or Mathematica.

@Carl Love Thanks a lot for that explanations on the second-derivative notation. It certainly helps a lot.

For the B[1,6] and B[2,6] part, I removed them because they are constants and they equal to zero

@tomleslie 

Thanks for your explanations.

Yes,  w(x,y)=0 but only at the edges of the rectangle, which explain the set of bc(s) but not in the inside of the rectangle

bc1 := w(0,y) = 0; # @ x=0 edge;
bc2 := w(a,y) = 0;  # @ x=a edge;
bc3 := w(x,0) = 0; # @ y=0 edge;
bc4 := w(x,b) = 0; # @ y=b edge;

The objective is to obtain w(x,y) in the in the inside of the rectangle which is not zero

@Carl Love 

Once again thanks for taking your time to help out. 

Based on your initial comment, I removed the constant terms since they play no role, which simplify to D[2](w)(x,b) = 0. I also change two of those D[2]s to be D[1]. but I was still getting the error message.

I included the original boundary conditions may be that way you will be able advise me on a better way to go about it 

@Carl Love

Thanks Dr. Carl. I tried that as well but I'm still getting the error message. 

@tomleslie 

Exactly, w(x, y) = 0 for pdsolve(sol, w(x,y)); which is not what I want. However, when I used pdsolve(myPDE1), it gives me w(x,y) with some constant terms. I notice for example bc1 := w(0, y) = 0 is like a boundary condition at a point whereas, it was suppose to be for the edge. Is there a way to specify the boundary conditions for the edges of the rectagular domain?