Download mw.mw

restart;
f:=n->piecewise(n=0,1,n>=1,sum('f(k)',k=0..n-1));
seq(f(i), i=0..10);
        1, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512

You have

plot(l, ....);

What is l? It is not defined as far as I can tell.  Perhaps you wanted to say

plot(x[1], ....);

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.