I may not be understanding your question correctly, but perhaps the following can help.

Download FourierSeries.mw

Solving this PDE is a textbook example an application of the Fourier series.  Tomleslie has already sketched the first half of the process.  In the attached worksheet I have carried it out all the way to the end.  There is one correction, however.  Tomleslie's calculations end up with X(x) expressed as an exponential function and Y(y) as a trigonometric function.  That won't work.  It should be the other way around in order to satisfy the boundary conditions bc[3] and bc[4].  The solution turns out to be

which agrees with what Maple 2017's pdsolve() produces, as shown by Kitonum.

Here's what the solution looks like:

Worksheet:   mw.mw

Update: The previous worksheet had a couple of careless typos, although the final solution was correct.  In this new version I have fixed the typos.

Worksheet ver 2: mw-ver2.mw

@Bendesarts After fixing the brackets issue we can make some progress but the results are not very encouraging.  Let's assume all variables other than g are positive:

indets(TableauRouth) minus {g}:
(x -> x>0)~(%):
rels := %[];
                  rels := 0 < ca, 0 < ka, 0 < keqc, 0 < ma, 0 < meqc, 0 < p

and then try the easier task of solving the equality TableauRouth[4,1]=0 for g:

solve(TableauRouth[4,1]=0, g) assuming rels;

If you try that, you will see that Maple does succeed in solving the equation but the result is so huge as to being useless. This tells us that there is no simple necessary and sufficient condition in terms of g for the stability of your system.

You should refer to matrix entries through square brackets, not parentheses.

Change TableauRouth(2,1)  to TableauRouth[2,1] to fix the Maple coding problem.  After that you will have to deal with the mathematics of your question, and that may not be very easy.

Here is how:

See Maple's help on fieldplot3d.  This is an example from there:

with(plots):
fieldplot3d([2*x,2*y,1],x=-1..1,y=-1..1,z=-1..1,grid=[5,5,5]);

Your argument here is backwards: You should check the correctness of Matlab with the help of Maple, not the other way around.

Matlab does its calculations in the machine hardware which carries about 16 digits of precision.  Maple can do its floating point calculations in software with any degree of precision.  The default is 10 digits which is less than Matlab's, but you may request higher.  For instance:

K := 25888158.4218766;
K-sqrt(K^2-1);              # default 10 digits precision
    0.
evalf(K-sqrt(K^2-1), 30);   # 30 digits precision
    1.93138496702600*10^(-8)

Compare this correct answer to Matlab's not-so-good answer.

You may change Maple's default precision of 10 by setting the value of Digits , as in:
Digits := 30;
Then all floating point calculations will be done with 30 digists of accuracy in software.  This will slow things down, therefore I don't recommend doing that unless you really have a need for it.

As an alternative to what kitonum has suggested, you may try:

ss1:=-2:  ss2:=1/2:  ss3:=-5:
printf("The 3 equation in A, B, C for s=%a, s=%a, and s=%a:\n", ss1, ss2, ss3);
     The 3 equation in A, B, C for s=-2, s=1/2, and s=-5:

Here is an alternative to what rlopez has suggested.

restart;
with(Physics):
Typesetting:-Settings(typesetdot=true):

Here I will pick a Lagrangian function just for illustration.  Replace it with yours.

L := 1/2*a[1]*diff(q[1](t),t)^2 + 1/2*a[2]*diff(q[2](t),t)^2
   + b*diff(q[1](t),t)*diff(q[2](t),t)
   + c[1]*q[1](t) + c[2]*q[2](t);

The two Euler-Lagrange equations are:

de[1] := diff(diff(L, diff(q[1](t),t)), t) - diff(L, q[1](t)) = Q[1](t);
de[2] := diff(diff(L, diff(q[2](t),t)), t) - diff(L, q[2](t)) = Q[2](t);

Here is an alternative to the solution presented by vv and kitonum.  It involves a little bit of a "cheat" in embedding the drawings in 3D.  The method can handle non-convex domains with equal ease.  Here are a couple of samples extracted from the worksheet mw.mw.

You are attempting to solve your equations in two phases: first for p and q, and then for M.  It is easier to solve all of them together at once.

There are four unknowns, p1, p2, q1, q2, from the first set of equations.  There are 16 knowns, M[i,j], from the second set of equations.  Put them together and solve a system of 20 equations in 20 unknowns, and then you are done.

You didn't say what your differential equations are, therefore I will make up some arbitrary ones here.
de1 := diff(x(t),t,t) = -x(t) - diff(x(t),t) + y(t);
de2 := diff(y(t),t,t) = -y(t) - diff(y(t),t) + x(t);
ic := x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=1;

As you see, I have isolated the second derivative terms on the left-hand sides of de1 and de2.

Now solve the system:
dsol := dsolve({de1, de2, ic}, numeric);

Let's say we want to examine the solution at time t=1.23.  We do
dsol(1.23);

The accelerations at time t=1.23 are obtained from
eval(rhs(de1), dsol(1.23));
eval(rhs(de2), dsol(1.23));

This is only a slight variation of a question that you asked a month ago.  Here is how to solve it:

Worksheet: mw.mw