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}

If you want to use your rule to specific real numbers, then  write  n::realcons .

In general case write  n::anything .

Examples:

applyrule( n::realcons*F(K,L) = F(n*K,n*L), Pi*F(K,L) );

applyrule( n::anything*F(K,L) = F(n*K,n*L), k*F(K,L) );

                                   F(Pi*K, Pi*L)

                                    F(k*K, k*L)

Example 1. 

Find the least positive root of the equation  tan(x)+a*x=0  for  a  in the list  [0.1, 0.2, ...,0.9, 1]  (total 10 numbers).

restart;

map(fsolve, [seq(tan(x)+a*x=0, a=[seq(0.1*n, n=1..10)])], x=Pi/2..3*Pi/2);

assign(seq(x||i=%[i], i=1..10));

x1, x3;

   [2.862772588, 2.653662400, 2.498399006, 2.380644485, 2.288929728, 2.215707303, 2.155983633, 2.106375674, 2.064528178, 2.028757838]

                                                     2.862772588, 2.498399006

Example 2.

 Find the first  10  positive root of the equation  tan(x)+0.1*x=0 . As the parameter we mean the range in which we look for the root.

restart;

map2(fsolve, tan(x)+0.1*x=0, [seq(x=-Pi/2+Pi*n..Pi/2+Pi*n, n=1..10)]);

assign(seq(x||i=%[i], i=1..10));

x2, x5;

[2.862772588, 5.760557933, 8.708313831, 11.70267808, 14.73347234, 17.79083538, 20.86723817, 23.95736763, 27.05755029, 30.16523655]

                                                      5.760557933, 14.73347234

We have used the obvious fact that in each range  x=-Pi/2+Pi*n .. Pi/2+Pi*n  (n is integer)  the equation  tan(x)+a*x=0  (a>0)    has a unique root.

Default Maple operates with angles expressed in radians:

convert((58+1/60*28+1/60^2*18)*degrees, radians);

                                 (35083/108000)*Pi

Sys1 := {diff(r(t),t) = r(t)^2*sin(theta(t)), diff(theta(t),t) = -r(t)^2*(-2*cos(theta(t))^2+1)}

DEtools[phaseportrait](Sys1, [r(t),theta(t)],t=-1..1,[[r(0)=1,theta(0)=0]]);