Here is my code

restart;

with(LinearAlgebra):

P := unapply(Matrix(4, 4, {(1, 1) = p[w], (1, 2) = p[x], (1, 3) = p[y], (1, 4) = p[z], (2, 1) = -p[x], (2, 2) = p[w], (2, 3) = -p[z], (2, 4) = p[y], (3, 1) = -p[y], (3, 2) = p[z], (3, 3) = p[w], (3, 4) = -p[x], (4, 1) = -p[z], (4, 2) = -p[y], (4, 3) = p[x], (4, 4) = p[w]}),(p[w],p[x],p[y],p[z]));

evalm(P(x,y,z,w));

(y1, y2, y3, y4) -> rtable(1 .. 4, 1 .. 4, {(1, 1) = y1,

(1, 2) = y2, (1, 3) = y3, (1, 4) = y4, (2, 1) = -y2,

(2, 2) = y1, (2, 3) = -y4, (2, 4) = y3, (3, 1) = -y3,

(3, 2) = y4, (3, 3) = y1, (3, 4) = -y2, (4, 1) = -y4,

(4, 2) = -y3, (4, 3) = y2, (4, 4) = y1}, datatype = anything,

subtype = Matrix, storage = rectangular, order = Fortran_order)

S:=Matrix(3,5):

for j from 1 to 3 do;

for k from 1 to 5 do;

S[j,k]:=P(seq(RandomTools[Generate](integer(range = 0 .. 4)),i=1..4));;

end do:

end do:

S;

W:=Matrix(3,5):

for j from 1 to 3 do;

for k from 1 to 5 do;

W[j,k]:=Transpose(S[j,k]);;

end do:

end do:

W;

M:=Matrix(3,3,shape=symmetric):

for j from 1 to 3 do;

for k from 1 to j do;

M[j,k]:=P(x[j,k],y[j,k],z[j,k],w[j,k]);

end do:

end do:

M;

N:=Matrix(3,3,shape=antisymmetric):

for j from 1 to 3 do;

for k from j+1 to 3 do;

N[j,k]:=P(u[j,k],v[j,k],q[j,k],r[j,k]);

end do:

end do:

N;

H:=M+N;

(Transpose(W)[2].H.Transpose(Transpose(S)[2]));

I would like the output to be a 4x4 matrix.