Perhaps this?

> allfactors:= P -> 
     map(convert, combinat[choose](map(t -> t[1]$t[2],factors(P)[2])),`*`);

Is this something like what you want?

> map(factor, collect(%, [eta, g], distributed));

(1/9)*(nu-1)^2*(nu-2)*(nu-3)*g^2*eta^4+(1/3)*(2*(nu-1))*g^2*eta^2+(1/6)*(nu-1)*(lambda*nu^2-2*nu*lambda-2-3*lambda)*eta^2

> P:= plot(...):
  plots[display](P, labels = [AAAA, BBBB], 
     labeldirections= [horizontal, vertical]);

It is very likely that the antiderivative can't be expressed in "closed form", and there doesn't appear to be anything special about the endpoint 0.2 that might lead one to think that special techniques (such as residues) could be applied to the definite integral.  Perhaps you could use power-series solutions, e.g. (with J your integral)

> sort(mtaylor(J, [a,b,c,d]), [a,b,c,d], ascending);

.2000000000+.1000000000*d^2+.4000000000e-1*c*d+.5333333333e-2*c^2+.8000000000e-2*b*d+.2400000000e-2*b*c+.2880000000e-3*b^2+.1600000000e-2*a*d+.5120000000e-3*a*c+.1280000000e-3*a*b+.1462857143e-4*a^2-.2500000000e-1*d^4-.2000000000e-1*c*d^3-.7999999999e-2*c^2*d^2-.1600000000e-2*c^3*d-.1280000000e-3*c^4-.4000000000e-2*b*d^3-.3600000000e-2*b*c*d^2-.1152000000e-2*b*c^2*d-.1280000000e-3*b*c^3-.4320000000e-3*b^2*d^2-.2880000000e-3*b^2*c*d-.4937142857e-4*b^2*c^2-.2468571428e-4*b^3*d-.8640000000e-5*b^3*c-.5759999999e-6*b^4-.8000000000e-3*a*d^3-.7680000000e-3*a*c*d^2-.2560000001e-3*a*c^2*d-.2925714286e-4*a*c^3-.1920000000e-3*a*b*d^2-.1316571428e-3*a*b*c*d-.2304000000e-4*a*b*c^2-.1728000000e-4*a*b^2*d-.6144000000e-5*a*b^2*c-.5529600000e-6*a*b^3-.2194285714e-4*a^2*d^2-.1536000000e-4*a^2*c*d-.2730666666e-5*a^2*c^2-.4096000001e-5*a^2*b*d-.1474560000e-5*a^2*b*c-.2010763636e-6*a^2*b^2-.3276800000e-6*a^3*d-.1191563636e-6*a^3*c-.3276800000e-7*a^3*b-.2016492308e-8*a^4

Hmm, interesting: the <maple> tag doesn't preserve the order of terms in the sum.  That should have been

.2000000000+.1000000000*d^2+.4000000000e-1*c*d+.5333333333e-2*c^2+.8000000000e-2*b*d+.2400000000e-2*b*c+.2880000000e-3*b^2+.1600000000e-2*a*d+.5120000000e-3*a*c+.1280000000e-3*a*b+.1462857143e-4*a^2-.2500000000e-1*d^4-.2000000000e-1*c*d^3-.7999999999e-2*c^2*d^2-.1600000000e-2*c^3*d-.1280000000e-3*c^4-.4000000000e-2*b*d^3-.3600000000e-2*b*c*d^2-.1152000000e-2*b*c^2*d-.1280000000e-3*b*c^3-.4320000000e-3*b^2*d^2-.2880000000e-3*b^2*c*d-.4937142857e-4*b^2*c^2-.2468571428e-4*b^3*d-.8640000000e-5*b^3*c-.5759999999e-6*b^4-.8000000000e-3*a*d^3-.7680000000e-3*a*c*d^2-.2560000001e-3*a*c^2*d-.2925714286e-4*a*c^3-.1920000000e-3*a*b*d^2-.1316571428e-3*a*b*c*d-.2304000000e-4*a*b*c^2-.1728000000e-4*a*b^2*d-.6144000000e-5*a*b^2*c-.5529600000e-6*a*b^3-.2194285714e-4*a^2*d^2-.1536000000e-4*a^2*c*d-.2730666666e-5*a^2*c^2-.4096000001e-5*a^2*b*d-.1474560000e-5*a^2*b*c-.2010763636e-6*a^2*b^2-.3276800000e-6*a^3*d-.1191563636e-6*a^3*c-.3276800000e-7*a^3*b-.2016492308e-8*a^4

The question must be restated as: are the three components of the difference of the two vectors ever all zero?  
You could try something like this:

> fsolve(add((r(t)[i]-s(t)[i])^2, i=1..3), t = - 2 .. 2);

If the result is a number, or a sequence of numbers, the answer is yes.  If there is no result or an unevaluated call to fsolve, the answer is no.

1) [...] in Maple makes a list, or an index for a table, list etc.  For grouping, use (...).
2) The complex number sqrt(-1) is I, not i.
3) You need multiplication signs * (or space in 2D Math input) between the parenthesized expressions in M1.

With those corrected:

> simplify(M1);

-M*(b+sin(theta)*a*exp(theta*I)*I-cos(theta)*a*exp(-I*theta))/(-b+a)

> simplify(convert(%,trig));

-M/(-b+a)*(2*I*sin(theta)*a*cos(theta)+b-a)

 > combine(%,trig);

(-M*b-I*M*a*sin(2*theta)+M*a)/(-b+a)

You might try using piecewise instead of `if` in your differential equation.  But it really makes no sense to have a quantity in a differential equation have special values at only a finite number of points, different from what it is everywhere else.   And it especially makes no sense for the value everywhere else to be NULL.  I really don't understand what you're trying to accomplish with this.

1) Check your spelling. 

2) C2 is not the same as c2.

3) You should use Matrix rather than matrix.

4) You should assign a value to finalmatrix[n,m] after your loop on q, instead of assigning a value to finalmatrix during that loop.

5) You want the procedure to return finalmatrix.

6) "The dimensions n, m, q should be
input parameters in your procedure."

For example:
 

> with(Statistics):
  X:= RandomVariable(Normal(0,1));
  f:= PDF(X,x);
  x1:= Quantile(X,0.05, numeric);
  x2:= Quantile(X,0.95, numeric);
  plots:-display([plot(f,x=-3..3),
                  plot(f,x=-3..x1,colour=gray,filled=true),
                  plot(f,x=x2..3,colour=gray,filled=true)]);

> seq(nops(L||i), i=1..3);

I'm assuming L1, L2, L3 are global variables.   This won't work for local variables or formal parameters of a procedure.

The equations (or rather expressions: I assume the equations are EQ1 = 0, ..., EQ4 = 0) you gave us don't have any a2 in them: the variables are a, P, A, C2 and alfa.  Moreover, they are linear in P, A, C2 and alfa. 

> solve({EQ1,EQ2,EQ3,EQ4},{P,A,C2,alfa});

{A = .2000000000e-1*(.5968744856e29*a^3+.4051173415e22-.6803542549e19*a^2)/(.1378851290e18*a^2+.8210378122e20+.6729306195e19*a), C2 = -.4000000000e-13*(-.5011956467e38*a^3+.6803542550e31-.1834515619e40*a^2+.5576250004e30*a+.1026962134e37*a^4)/(.1378851290e18*a^2+.8210378122e20+.6729306195e19*a), P = .4000000000e-9*(.2050351869e28*a^3+.2783275608e21+.9348468735e18*a^2+.2281199904e20*a)/(.1378851290e18*a^2+.8210378122e20+.6729306195e19*a), alfa = .4000000000e-14*(.6803542549e32*a^2+.2025586708e35+.1660190673e34*a+.8558017830e36*a^6-.7643815078e40*a^4-.1252989117e39*a^5+.1858750001e31*a^3)/(.1378851290e18*a^2+.8210378122e20+.6729306195e19*a)}

I am quite confident that this is correct (up to roundoff error).  Yes, there will be some large values for the variables, e.g. for a = 40 you get

{A = .1335914871e12, C2 = 245770363.6, P = 91.78129128, alfa = -.2020898993e12}

But if these results are very different from the "real ones", it can only be because the equations are wrong.

Well, you don't need all those x[i] <= x clauses:
piecewise(x < x[1], hz[1], x < x[2], hz[2], etc)
would do.  But if you want to use x as a symbolic variable, it's probably better to use a different letter for the indexed object, maybe X[i] rather than x[i]. 

And for producing it, you could use something like

> piecewise(op(ListTools[Flatten]([seq]([x < X[i], hz[i]],i=1..10)))); 

Suppose the probability density for trait A in group B (i.e. the conditional density given that the individual is in group B) is f_B(x),  the probability density in the complement of group B is f_C(x), and the probability of being in group B is p, Then the overall probability density for trait A is f_all(x) = p*f_B(x) + (1-p)*f_C(x).  Therefore f_C(x) = (f_all(x) - p*f_B(x))/(1-p). 

Caution: if two of these densities are normal, the third will _not_ be normal, except in trivial cases.

The bug is still present in 13.02.

DEplot uses the differential equation, not the solution.  If you're using dsolve to get the solution, you just need plot to plot the result.  Thus:

> plot(rhs(eq2), t = 0 .. 25);

You could use DEplot this way:

> with(DEtools):
  DEplot(subs(lambda=.1, eq1), n(t), t = 0 .. 25, [n(0)=10],
    arrows = none);