Equality of matrices and arrays is tested by the command  LinearAlgebra[Equal] :

c := Array([1, 2, 3]):

d := Array([1, 2, 3]):

LinearAlgebra[Equal](c, d);

                true

Such elements may be multiple.

Example:

L:=[1.2, 0.25, 5.6, 7.2, 0.5, -0.25, -4, 6]:

L1:=map(abs, L);

a:=min(L1);

ListTools[SearchAll](a, L1);

If  M  is the number of the rows, the last number in the first column   -1/(2*M*(M-2)  is not in accord with the first numbers in the odd-numbered rows. So I think that the number  2  in the denominator is superfluous. If this 2 removed, then we obtain the solution:

S:=proc(k::{symbol,realcons}, N::{odd,positive})

``(sqrt(2)/2^k)*``(Matrix(N, {seq((m,1)=`if`(m::even,0,-1/m/(m-2)), m=1..N)}));

end proc:

Example:

S(k, 7);

With Maple it will be  shorter to use  maximize  command:

restart:

h:= r-> 2*sqrt(100-r^2):

V:= r-> Pi*r^2*h(r):

maximize(V(r), r=0..10, location);

assign(op(%[2])[1]);

radius=r, height=h(r); 

As regards the notion of equivalence of two matrices, then look here . 

Two matrices are equivalent if and only if they have the same size and the same rank. So the code be very simple:

is(LinearAlgebra[Rank](A) = LinearAlgebra[Rank](B) and op(A)[1 .. 2] = op(B)[1 .. 2]);

Edited.

I did not find errors in the generated matrix. This matrix depends on the list  w . Maybe you mean the following order of unknowns:

w := [seq(seq(u[i,j],j=1..N),i=1..N)];

w := [u[1, 1], u[1, 2], u[1, 3], u[1, 4], u[2, 1], u[2, 2], u[2, 3], u[2, 4], u[3, 1], u[3, 2], u[3, 3], u[3, 4], u[4, 1], u[4, 2], u[4, 3], u[4, 4]]

Maybe you mean the diagonal matrix of n blocks?

Example:

A:=Matrix(3, [1,2,3,2,1,3,3,2,1]);

LinearAlgebra[DiagonalMatrix]([A $ 3]);

remove(has, [x+p*y, x+y, y+z], p);

                     [x+y, y+z]

restart;

solve(eval({z = x^2+2*y^2, z = -2*x^2-y^2+3, -2*x^2-y^2+3 = x^2+2*y^2}, x = t), {y, z});

allvalues(%);

assign(%[1], x = t);

4*(int(sqrt((diff(x, t))^2+(diff(y, t))^2+(diff(z, t))^2), t = 0 .. 1));

evalf(%);

Suppose we know the number of elements of a set  A  in which the numbers of elements  its  N  subsets and their intersections in two, in three and so on are known. The procedure  None  computes the number of elements  of  A belonging none of these subsets. 

None:=proc(L::list(posint), N::posint)

local K;

uses combinat;

K:=choose(N);

L[1]-add(L[k]*(-1)^(nops(K[k])-1), k= 2..nops(L));

end proc:

Example:

None([800, 224, 240, 336, 64, 80, 40, 24], 3);

                                 160

800 - (224 + 240 + 336 - (64 + 80 + 40) + 24);

                                    160

At the point  x = 0, violated the conditions of existence and uniqueness of solutions. You can find a unique solution if will take a close point and will be solved numerically.

restart;

ode := (-6 + 3*x - 3*x^2 + 2*x^3)*y(x) + x*(6 - 3*x + x^3)*diff(y(x),x) + x^2*(-3 + 3*x - 3*x^2 + x^3)*diff(y(x),x$2)= 0;

ic := y(0.001)= 0, D(y)(0.001)= 1:

sol:=dsolve({ode, ic}, numeric);

plots[odeplot](sol, [x, y(x)], 0.001 .. 1);

The task is easy to solve by using Banach fixed-point theorem (contraction mapping principle). See  here

Consider the function   f:=x->1/2*(x+2/x) . The segment  [1, 2]  is a complete metric space in standard metric. 

Enough to check two conditions:
1) The function  f  maps the segment  [1, 2]  into itself.
2) The mapping  f  is a contraction.

This is easily done manually or with Maple.

Checking with Maple:

f:=x->1/2*(x+2/x):

is(f(x)>=1 and f(x)<=2) assuming x>=1,x<=2;  # Checking of first condition

abs('f'(a)-'f'(b))<=maximize(abs(diff(f(x),x)), x=1..2)*abs(a-b);  # Checking of second condition

The limit of the sequence  {u(n+1)=1/2*(u(n)+2/u(n)), u(0)=1}  is the fixed point of mapping  f  ;

solve({x=f(x), x>=1, x<=2});

Procedure  P  solves your problem for any number of couples  (m is the number of couples).

N:=proc(m::posint)

local P, C1, C, P1;

P:=(2*m)!;  # Total number of variants

C1:=2*binomial(2*m-1,1)*(2*m-2)!;  # Number of variants in which the members of the first couple are close

C:=k->binomial(2*m-k,k)*2^k*k!*(2*m-2*k)!;  # C(k) is the number of variants in which the members of the first  k couples are close

P1:=binomial(m,1)*C1+add(binomial(m,k)*(-1)^(k-1)*C(k), k=2..m); # Number of variants in which the members at least one couple are close (the exclusion-inclusion principle)

P-P1;

end proc:

Examples:

seq(N(m), m=1..20);  # The last number is the answer to your problem

Addition: independent check by the brute force method for  m=3 :

N:=0:  # {a,b}, {c,d}, {e,f} are three couples

for i in combinat[permute]([a,b,c,d,e,f]) do

if {i[1],i[2]}<>{a,b} and {i[1],i[2]}<>{c,d} and {i[1],i[2]}<>{e,f} and

{i[2],i[3]}<>{a,b} and {i[2],i[3]}<>{c,d} and {i[2],i[3]}<>{e,f} and

{i[3],i[4]}<>{a,b} and {i[3],i[4]}<>{c,d} and {i[3],i[4]}<>{e,f} and

{i[4],i[5]}<>{a,b} and {i[4],i[5]}<>{c,d} and {i[4],i[5]}<>{e,f} and

{i[5],i[6]}<>{a,b} and {i[5],i[6]}<>{c,d} and {i[5],i[6]}<>{e,f} then N:=N+1 fi:

od:

N;

                                                             240

assume(a::realcons, b::realcons):

eq:=I*(a+3*b)+2*b+5*a = 3+2*I:

solve({Re(lhs(eq))=Re(rhs(eq)), Im(lhs(eq))=Im(rhs(eq))});

                              {a = 5/13, b = 7/13}