The following method works:

eq2:=subsop(seq([1,k]=0, k=2..4), [1,5]=q11, eq1);

If you want to build a pyramid of arbitrary, it is convenient to use a special procedure . Formal arguments are: X - a list of names for the vertices of the pyramid, Y - a list of the coordinates. After specifying the pyramid for it's building  plots[display] command should be used. You can add other objects.

Code of Procedure:

Pyramid := proc(X::(list(symbol)), Y::(list(list)))

local L, L1, T, T1, n;

n:=nops(X);

L:=plottools[polygon]([seq(Y[i],  i=1..n-1)], style=line, color=blue, thickness=2);

L1:=seq(plottools[line](Y[i],Y[n], color=blue, thickness=2), i=1..n-1);

T:=plots[textplot3d]([seq([op(Y[i]), X[i]], i=1..n-1)], align={left,below}, color=black, font=[TIMES,ROMAN,18]);

T1:=plots[textplot3d]([op(Y[n]),X[n]], align=above, color=black, font=[TIMES,ROMAN,18]);

L, L1, T, T1;

end proc:

Example of constructing a pyramid with a base:

plots[display](Pyramid([A,B,C,D,S],[[0,0,0],[2,0,0],[2,3,0],[0,3,0],[0,0,4]]),plot3d(0,x=-1..3,y=-1..4,style=surface,color=pink), scaling=constrained, orientation=[35, 65]);

Use subs command instead of algsubs command.

print(`Answer:`); print(sqrt(5));

I do not know how to write it in a single line.

L:=combinat[permute]([0$8, 1$8], 8):

V:=[seq(v||i, i=1..8)]:

M:=[seq(add(L[k,i]*V[i], i=1..8), k=1..2^8)];

M is the list of all required sums.

This problem "to find the minimal set of the edges of G, such that its annihilation transforms G to a planar graph" is known as Graph Planarization Problem. It is proved that this problem is NP - complete. Therefore, there is no exact  algorithms for solving for large graphs,  only the approximate algorithms.

The presented procedure MAX is not exact.

The example. Consider the graph

G:=Graph(7, {{1, 2}, {1, 4}, {1, 5}, {1, 6}, {1, 7}, {2, 4}, {2, 6}, {2, 7}, {3, 4}, {3, 5}, {3, 6}, {3, 7}, {4, 5}, {4, 6}, {4, 7}, {5, 7}, {6, 7}});

Applying the procedure MAX, we obtain a planar graph with 14 edges, ie 3 edges removed. In fact, enough to remove the two edges {1, 2} and {1, 6} .

I ran your code, but I could not wait until the end of calculations. The rewritten code works fine:

J := u -> evalf(Int(1/(sqrt(x^2+u^2)*(exp(sqrt(x^2+u^2))+1)), x = 0 .. infinity)):

plot([seq([Delta, exp(-2*J(Delta))], Delta = 0 .. 2.5, 0.1)], thickness=2);

Your formal error is that when you write  Matrix([[M1],[M2],[M3]]) , the Maple sees it as a column matrix with unknown numbers M1, M2, M3, and not as a 3 by 3 matrix. Unknown matrix M of size 3 by 3 can be written so

M:=Matrix(3, 3, symbol=m);

and then you may refer to its elements in a known manner.

series(cos(x^2+x), x=0, 11):

a:=convert(%, polynom):

b:=[seq(coeff(a, x, n), n=0..10)];

The procedure P solves your problem:

P:=proc(k, i)

local L, s, M, K;

if not (k>=2*i and 2*i>=4) then

error `should be k>=2*i>=4` fi;

if not type(k, posint) or not type(i, posint) then

error `should be k is integer and i is integer` fi;

L:=[];

for s from k-2*i to k-i do

M:=combinat[composition](s+i, i);

K:=seq([seq(M[j,l]-1, l=1..i)], j=1..nops(M));

L:=[op(L), K];

od;

L;

end proc;

Example:

P(8, 3);

Matrix А isn't singular:

LinearAlgebra[Determinant](A);

                 0.03893229110      

with(geom3d):

Cd:=coordinates(M):

subs(x=Cd[1], y=Cd[2], z=Cd[3], Equation(P));

Объём_части_эллипсои.mws

Read more on the Russian site  http://forum.exponenta.ru/viewtopic.php?t=12179  and in the attached file 

For the functional coloring all surfaces should be given by explicit equations, rather than using graphical primitives.

A := plot3d([u, v, 1], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2+1)):

B := plot3d([u, 0, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

C := plot3d([0, u, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

E := plot3d([1, u, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

F := plot3d([u, 1, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2+1)):

G := plot3d([u, v, 0], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

plots[display](A, B, C, E, F, G, view = [-.27 .. 1.27, -.27 .. 1.27, -.27 .. 1.27], lightmodel = light4, orientation = [20, 70]);

It is clear that points equidistant from the origin have the same color. Detailization of coloring depends on numpoints option.

plot3d([ sqrt(5+4*sin(t))*sin(s), sqrt(5+4*sin(t))*cos(s),cos(t)], t=0..2*Pi, s=0..2*Pi, scaling=constrained, color=(5+4*sin(t))^2+cos(t)^2, axes=normal, style=surface,  view=[-3.7..3.7, -3.7..3.7, -1.7..1.7], numpoints=10000, lightmodel=light4, orientation=[20, 60]);