The penultimate line of code  

M[1..6,1..6]:=Matrix(3,3,shape=identity)

is not evaluated because it is written in text mode.

Yes, it's possible. Use the formula of differentiation of parametric function  Diff(y(x), z)=Diff(y(x), x)/Diff(z(x), x)

Example:

y:=x->sin(x): z:=x->x-y(x):

diff(y(x), x)/diff(z(x), x);

f:=x->piecewise(x>=2 and x<=6, sqrt(x), undefined);

D(f)(1);

plot(f, 0..7, 0..3, scaling=constrained);

restart;

Dist := proc (A::list, B::list)

sqrt((A[1]-B[1])^2+(A[2]-B[2])^2)

end proc:  # Procedure for calculation of the distance between two points defined by lists

A := [1, 1]: B := [x1, 2]: C := [x2, y2]:  # The centers of the circles

solve({Dist(A, B) = 3, Dist(A, C) = 1+y2, Dist(B, C) = 2+y2}, useassumptions) assuming x1 > 0, x2 > 0, y2 > 0;

assign(%);

Circle1 := (x-1)^2+(y-1)^2 = 1;

Circle2 := (x-x1)^2+(y-2)^2 = 4;

Circle3 := (x-x2)^2+(y-y2)^2 = expand(y2^2);

plots[implicitplot]([Circle1, Circle2, Circle3], x = -1 .. 6, y = -1 .. 5, thickness = 2, color = [blue, green, red], scaling = constrained, numpoints = 10000);

For older versions (at least for <=12)  ~  command doesn't work. The adjustment will be

seq( x^3-3*x^2-24*x+8, x=[-2,4,6,9]);

                             36, -72, -28, 278

In new versions this variant works also.

M:=3:  N:=3:

H:=1/N*Matrix(M,{seq((i,1)=t[f]^i/i, i=1..M)}):

E:=1/N*Matrix(M, {seq((i,i+1)=1/i, i=1..M-1)}):

P:=<seq(`<|>`(Matrix(M) $ M-1-k,E,H $ k),k=M-1..0,-1)>;

restart;

N := 8:

a := -Pi:

b := 3 * Pi:

data := []:

for n from 0 to N do

x[n] := a + n * (b - a)/N;

data := [op(data), [x[n], sin(x[n])]];

od:

data;

plot(data);

If the list is long better to write the following equivalent but more effective variant:

restart;

N := 1000:

a := -Pi:

b := 3 * Pi:

for n from 0 to N do

x[n] := a + n * (b - a)/N;

data[n] := [x[n], sin(x[n])];

od:

data := convert(data, list):

plot(data);

 Addition:  Instead a loop  shorter to use  seq  command:

restart;

a:=-Pi: b:=3*Pi:

N:=1000:

data := [seq([k, sin(k)], k=[a + n * (b - a)/N $ n=0..N])]:

plot(data);

Because we already know that  sum(1/n^2, n=1..infinity)=Pi^2/6 , the same can be written slightly shorter:

floor(fsolve(sum(1/n^2, n= N..infinity)=.001));

                              1000

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];