It looks better if only one parameter is changed.

For thickness (every value of the thickness takes about 2 seconds):

plots[display]([seq(plot(exp(-x^2)*sin(Pi*x), x = -2 .. 2, color = red, thickness = t, caption = typeset("Thickness =  ", t), titlefont = [TIMES, BOLD, 16]) $ 20, t = 1 .. 15)], insequence = true);

The same for colors (7 colors of the rainbow):

plots[display]([seq(plot(exp(-x^2)*sin(Pi*x), x = -2 .. 2, color = c, thickness = 4, caption = typeset("Color =  ", c), titlefont = [TIMES, BOLD, 16]) $ 20, c = [red, orange, yellow, green, blue, "Indigo", violet])], insequence = true);

Of course, these simple examples easier to solve by hand. 

Solution with Maple:

solve(3*x-5>=0);

solve(x+2<0 or x+2>0);

1)  nops(T)  is the number of elements of the list  T . For your example  

nops(T)=binomial(512, 3)=512*511*510/3!=22238720

2) Your code is very inefficient due to line  T := [op(T), [i,j,k]];

The more effective  and short way to use the code:

T:=combinat[choose]([$ 1..512], 3):

3) I do not recommend you to display this list, because it is very long. Can display its individual elements;

T[1],  T[1000000],  T[10000000];

        [1, 2, 3], [8, 274, 499], [93, 148, 227]

v:=<1,-2,3,-4>;  n:=op(1,v);

Matrix(`<|>`(v $ n));

1) I don't know what is gap system.

2) In Maple:

[op(x+2*y+z^2)];

       [x, 2*y, z^2]

1) It's even more simple task. First, we create the list  S  to iterate through, then  we pass along the list in a single loop:

restart;

n:=0:

S:=combinat[permute]([$ 3..14], 5):

for s in S do

F, H, J, K, L:=op(s):

if 1+F-J > 0 and 15+H-K > 0 and 18-J-K > 0 and 33-K-L > 0 and 34-H-J-L > 0 and -15-F+J+K+L > 0 and 16-F-H+J+K > 0 then n:=n+1: M[n]:=[F,H,J,K,L]: fi:

od:

M:=convert(M,list):

nops(M);

M[1..30];

M[-30..-1];

2) DirectSearch  package requires you to download from  Maple Application Center

 http://www.maplesoft.com/applications/view.aspx?SID=101333

Of your conditions is easy to get bounds for all variables:  1<=J<=17 ,  1<=K<=17 ,  1<=L<=32 ,  1<=H<=32 ,  1<=F<=36 . Next we find all the solutions by the brute force method. Displayed the total numbers of solutions, the first 30 and the last 30 solutions:

restart;

n:=0:

for F to 36 do

for H to 32 do

for J to 17 do

for K to 17 do

for L to 32 do

if 1+F-J > 0 and 15+H-K > 0 and 18-J-K > 0 and 33-K-L > 0 and 34-H-J-L > 0 and -15-F+J+K+L > 0 and 16-F-H+J+K > 0 then n:=n+1: M[n]:=[F,H,J,K,L]: fi:

od: od: od: od: od:

M:=convert(M,list):

nops(M);

M[1..30];

M[-30..-1];

For  alpha=-2:

eq:=diff(x(t),t)=exp(2*x(t))+alpha-2*x(t):

A:=DETools[dfieldplot](eval(eq, alpha=-2), x(t), t=-2..2, x=-2..2,arrows=SLIM, color=blue):

x1:=0.5730966103: x2:=-0.9207028302:

B:=plot([x1, x2], t=-2..2, color=red, thickness=2):

T:=plots[textplot]([[1,x1+0.1,x=x1], [1,x2+0.1,x=x2]], font=[TIMES,ROMAN,14]):

plots[display](A,B,T); 

Thus we see that  x1  is a point of unstable equilibrium, and  x2  is a point of stable equilibrium

1) You should write  exp(2*x(t))  instead of  e^{2*x(t)} .  e  is just a symbol in Maple.

2) Equilibrium points of an autonomous equation  dx/dt=f(alpha, x)  with parameter  alpha   are the solutions of the equation  f(alpha, x)=0   with respect to  x .  The overall picture can be seen from the graph of the implicit function:

plots[implicitplot](exp(2*x)+alpha-2*x=0, alpha=-6..6, x=-3..3, thickness=2, numpoints=3000, view=[-6..1,-3..3]);

It is easy to prove that the extreme right point has coordinates (-1, 0).

3)  In the point  alpha=-1  we have a saddle-node bifurcation, because if  alpha>-1  then no solutions, if  alpha=-1  then 1 solution (x=0), if  alpha<-1  then 2 solutions. These solution can be found numerically for every  alpha, for example:

eq:=exp(2*x)+alpha-2*x=0:

fsolve(eval(eq, alpha=-2), x=0..infinity);

fsolve(eval(eq, alpha=-2), x=-infinity..0);

                    0.5730966103

                   -0.9207028302

Your first equation  T1=0  is easily reduced to a cubic equation by change  x=r^2 . Therefore, the original equation T1 = 0 will have two real roots if and only if the corresponding cubic equation has only one real positive root. General cubic equation  a*x^3+b*x^2+c*x+d=0  has only one real positive root if and only if two conditions are true (see  http://en.wikipedia.org/wiki/Cubic_function ):

discriminant  Delta=18*a*b*c*d-4*b^3*d+b^2*c^2-4*a*c^3-27*a^2*d^2<0  and  a*d<0

All real roots, you can find exactly (symbolically) by the  RealDomain[solve]  command.

Example:

RealDomain[solve](2*r^6-6*r^4+5*r^2-3=0);

Before plotting all constants should be specified. Of cause, you can build several graphs at one plot.

Example:

bX:=-1:  aX:=5:

plot([seq(seq((X_0-aX)/a+(X_0-bX)/b, a=1..3), b=1..3)], X_0=bX..aX, thickness=2, color=[red,blue,brown, yellow,pink,green,gold,cyan,grey]);

Another example (a, b are constants, aX, bX are changed):

restart;

a:=-1:  b:=5:

assign(seq(A[n]=plot((X_0-n)/a+(X_0-(n+5))/b, X_0=n..n+5, thickness=2, color=[red,blue,brown, yellow,pink][n]), n=1..5)):

plots[display](seq(A[n], n=1..5));

Just check correctness of these rules for each example.

Solution of (a):

T:=(a,b,c,d)->Matrix([[a,a^2,a^3],[b,c,d]]):

M1:=T(a1+a2,b1+b2,c1+c2,d1+d2);

M2:=T(a1,b1,c1,d1)+T(a2,b2,c2,d2);

LinearAlgebra[Equal](M1,M2);

The plotting depends on the choice of contours and locations for text. If the number of contours is equal to  n, then I divide the range ​​of the function on  n  parts and choose the values ​​of the function at the midpoints of these parts.

restart;

f := (x, y)->x^2+y^2-5;

m, M := minimize(f(x, 0), x = 0 .. 10),  maximize(f(x, 0), x = 0 .. 10);

n, d := 5, (M-m)/n;

C := seq(m+(1/2)*d+d*(i-1), i = 1 .. n);

P := plots[contourplot](f, -10 .. 10, -10 .. 10, color=brown, contours = [C]):

T := plots[textplot]([seq([solve(f(x, x) = C[i])[1]$2, C[i]], i = 1 .. n)], font = [TIMES, ROMAN, 14], align = right):

plots[display](P, T);

 See  Statistics[Regression]  help page. 

I could not find the options for the exact solution. But the exact solution can be easily obtained by the simple formulas:

with(Student[LinearAlgebra]):

a:=<-2|3|2>;  b:=<7|-3|-4>;

Pr:=(a.b)/(Norm(b)^2)*b;   # exact projection

evalf(Pr);   # approximate projection

c:=a-Pr;  # exact orthogonal complement

evalf(c);   # approximate orthogonal complement

Norm(c);  # exact norm of orthogonal complement

evalf(%);   # approximate norm of orthogonal complement

 Addition:  See http://en.wikipedia.org/wiki/Vector_projection