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

Firstly, are created a list of all pairs of diagonals for such squares, and then in a loop filled with the remaining entries. The total number of such  Latin squares is equal to 960. For illustration  3 squares are displayed .

restart:

L:=combinat[permute](5):

ND:=0: Diag:=[]:

for i in L do

A:=convert(i, set) minus {i[3]}:

M:=combinat[permute](A):

for j in M do

if i[1]<>j[1] and i[2]<>j[2] and i[4]<>j[3] and i[5]<>j[4] and i[1]<>j[4] and i[2]<>j[3] and i[4]<>j[2] and i[5]<>j[1] then ND:=ND+1: Diag:=[op(Diag),[[seq(i[k], k=1..5)], [j[1], j[2], i[3], j[3], j[4]]]]: fi:

od: od:

ND: # The number of pairs of diagonals

Diag: # The list of pairs of diagonals

B:={1, 2, 3, 4, 5}:

N:=0: List:=[]:

for i from 1 to nops(Diag) do

for a in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,1,2], Diag[i,2,2]} do

for b in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,1,3]} do

for c in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,2,4], Diag[i,1,4]} do

for d in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,1,1], Diag[i,2,1]} do

for e in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,1,3]} do

for f in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,2,5], Diag[i,1,5]} do

for g in B minus {Diag[i,1,3], Diag[i,1,1], Diag[i,2,1]} do

for h in B minus {Diag[i,1,3], Diag[i,1,2], Diag[i,2,2]} do

for l in B minus {Diag[i,1,3], Diag[i,2,4], Diag[i,1,4]} do

for m in B minus {Diag[i,1,3], Diag[i,2,5], Diag[i,1,5]} do

for p in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,1,1], Diag[i,2,1]} do

for q in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,1,3]} do

for r in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,2,5], Diag[i,1,5]} do

for u in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,1,2], Diag[i,2,2]} do

for v in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,1,3]} do

for w in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,2,4], Diag[i,1,4]} do

if nops({Diag[i,1,1],a,b,c,Diag[i,2,5]})=5 and nops({d,Diag[i,1,2],e,Diag[i,2,4],f})=5 and nops({g,h,Diag[i,1,3],l,m})=5 and nops({p,Diag[i,2,2],q,Diag[i,1,4],r})=5 and nops({Diag[i,2,1],u,v,w,Diag[i,1,5]})=5 and nops({Diag[i,1,1],d,g,p,Diag[i,2,1]})=5 and nops({a,Diag[i,1,2],h,Diag[i,2,2],u})=5 and nops({b,e,Diag[i,1,3],q,v})=5 and nops({c,Diag[i,2,4],l,Diag[i,1,4],w})=5 and nops({Diag[i,2,5],f,m,r,Diag[i,1,5]})=5 then N:=N+1: List:=[op(List), Matrix([[Diag[i,1,1],a,b,c,Diag[i,2,5]], [d,Diag[i,1,2],e,Diag[i,2,4],f], [g,h,Diag[i,1,3],l,m], [p,Diag[i,2,2],q,Diag[i,1,4],r], [Diag[i,2,1],u,v,w,Diag[i,1,5]]])]: fi:

od: od: od: od: od: od: od: od: od: od: od: od: od: od: od: od: od:

N;  # N -  the total number of squares,  List - the list of all squares

List[1], List[100], List[960];

A possible way :

with(LinearAlgebra):

M:=RandomMatrix(6);

G:=DiagonalMatrix(RandomVector(6));

M - DiagonalMatrix([seq(M[i,i], i=1..6)])+DiagonalMatrix([a$6]);

M - DiagonalMatrix([seq(M[i,i], i=1..6)])+G;

a:=2: b:=4: c:=7: d:=2:

if b>3 then z:=a:  y:=c elif b<3 then z:=d:  y:=a end if:

z;  y;

2

7

Add  numpoints  option:

plot3d(sin(sqrt(x^2+y^2))/sqrt(x^2+y^2), x=-3..3, y=-3..3, axes=normal, numpoints=3000);

Use

Optimization[Minimize]  command.