@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

restart;
P := <1,1,1> + s*<2,-1,-1> + t*<-1,3,2>;
Q := <1, -3, 5>/sqrt(35);
t_range := 0..1;
s_range := 0..1;
center := eval(P, {t=`+`(op(t_range))/2, s=`+`(op(s_range))/2});
plots:-display([
  plot3d(P, t=t_range, s=s_range, color=blue, style=surface),
  plots:-arrow(center, Q, color=red)
]);

_U1 is a "dummy" variable of integration, because the solution of your differential equation is expressed as a convolution.  You may replace _U1 by any other symbol you like, for example:
restart;
de := diff(x(t),t) = a*x(t) + b(t)*u(t);
dsolve({de, x(0)=alpha}, x(t), method=laplace);
subs(_U1=s, %);