A := Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = -1, (2, 2) = 2});

B := Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = -2, (2, 2) = 3});

C := Matrix(2);

M:=4:

P:=Matrix([seq([C$k,A,B$M-1-k],k=0..M-1)]);

L := [sin(x), sin(x)-x+x^2, cos(x)-sin(x)]:

X := [seq(-Pi+(1/16)*Pi*k, k = 1 .. 31)]:

Y := map2(eval, L, map(t->x = t, X)):

P := [seq(zip(`[]`, X, [seq(Y[i, j], i = 1 .. 31)]), j = 1 .. 3)]:

A := plot(L, x = -Pi .. Pi, color = [red, blue, green], legend = L):

B := plot(P, x = -Pi .. Pi, color = [red, blue, green], style = point, symbol = solidcircle, symbolsize = 10):

plots[display](A, B, view = [-Pi .. Pi, -2 .. 4], scaling = constrained);

First you have to specify the values of the constants  _C1  and  _C2 .

Example:

y(t) = _C1*exp(-1.*t)*sin(.57736*t)+_C2*exp(-1.*t)*cos(.57736*t):

assign(%);

_C1, _C2 := 2, 1:

plot(y(t), t=Pi..2*Pi);

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