Verification that the algebra of 4x4 matrices is the tensor product of H with H (where H is the algebra of quaternions). 

> with(LinearAlgebra):
 

First some auxiliary 2x2 matrices. These will be used to construct 4x4 matrices: 

> J:=Matrix([[0,-1],[1,0]]);
JD:=Matrix([[-1,0],[0,1]]);
I2:=IdentityMatrix(2);
JJ:=Matrix([[0,1],[1,0]]);
Z2:=Matrix(2,2,0);
 

J := Matrix(%id = 146852896) 

JD := Matrix(%id = 147109568) 

I2 := Matrix(%id = 147333768) 

JJ := Matrix(%id = 147526936) 

Z2 := Matrix(%id = 147854152) 

The matrices I4, IL,JL,KL  given below generate one copy of the quaternions
and the matrices I4, IR,JR,KR generate another copy of the quaternions.
(Note that if one replaces IR by -IR, JR by -JR and KR by -KR then I4, IR,JR,KR will generate
an algebra isomorphic to the opposite algebra of H (actually the same as H but with a different basis.)  

> IL:=Matrix([[J,Z2],[Z2,J]]);
JL:=Matrix([[Z2,JD],[-JD,Z2]]);
KL:=Matrix([[Z2,-JJ],[JJ,Z2]]);
IR:=Matrix([[-J,Z2],[Z2,J]]);
JR:=Matrix([[Z2,I2],[-I2,Z2]]);
KR:=Matrix([[Z2,-J],[-J,Z2]]);
I4:=IdentityMatrix(4);
 

IL := Matrix(%id = 148292412) 

JL := Matrix(%id = 148714148) 

KL := Matrix(%id = 149213768) 

IR := Matrix(%id = 149413744) 

JR := Matrix(%id = 146828076) 

KR := Matrix(%id = 148706032) 

I4 := Matrix(%id = 148089776) 

Verification that the matrices I4, IL,JL,KL generate the quaternions and the matrices I4, IR,JR,KR generate another copy of the quaternions.
It suffices to show that IL^2 = -I4, JL^2 = -I4, KL^2 = -I4 and IL.JL.KL = -I4. And similarly for R replacing L. This is verified  below. 

> IL^2 +I4, JL^2 +I4, KL^2 +I4,IL.JL.KL + I4;
IR^2 +I4, JR^2 +I4, KR^2 +I4,IR.JR.KR + I4;  
 

Matrix(%id = 148449368), Matrix(%id = 148449368), Matrix(%id = 148449368), Matrix(%id = 148449368) 

Matrix(%id = 147707972), Matrix(%id = 147707972), Matrix(%id = 147707972), Matrix(%id = 147707972) 

We now verify that the following 16 matrices are linearly independent and hence form a basis for the 4x4 matrices. Since the 4x4 matrices has dimension 16 they must therefore be a basis for the algebra of 4x4 matrices. 

> B:=[I4,IL,JL,KL,IR,JR,KR,
IL.IR,IL.JR,IL.KR,
JL.IR,JL.JR,JL.KR,
KL.IR,KL.JR,KL.KR] :
 

> A:=simplify(add(a[i]*B[i],i=1..16));
 

A := Matrix(%id = 147576200)
A := Matrix(%id = 147576200)
A := Matrix(%id = 147576200)
A := Matrix(%id = 147576200)
A := Matrix(%id = 147576200)
A := Matrix(%id = 147576200) 

> solve({seq(seq(A[i,j]=0,i=1..4),j=1..4)});
 

{a[14] = 0, a[11] = 0, a[16] = 0, a[13] = 0, a[15] = 0, a[12] = 0, a[10] = 0, a[2] = 0, a[5] = 0, a[4] = 0, a[3] = 0, a[1] = 0, a[8] = 0, a[6] = 0, a[7] = 0, a[9] = 0}
{a[14] = 0, a[11] = 0, a[16] = 0, a[13] = 0, a[15] = 0, a[12] = 0, a[10] = 0, a[2] = 0, a[5] = 0, a[4] = 0, a[3] = 0, a[1] = 0, a[8] = 0, a[6] = 0, a[7] = 0, a[9] = 0} 

We can write an arbitrary 4x4 matrix C as a linear combination of the basis B as follows: 

> C:=matrix(4,4,(i,j)->c[i,j]);
 

C := Matrix(%id = 149245748) 

Note that the coefficient of I4 is the trace of the matrix C. And it is known for example that trace(AB) = trace(BA).  

> sol:=op(solve({seq(seq(A[i,j]=C[i,j],i=1..4),j=1..4)},[seq(a[i],i=1..16)])):
for i from 1 to 16 do
print(sol[i], "=coefficient of", B[i]);
od;
 

a[1] = 1/4*c[1, 1]+1/4*c[4, 4]+1/4*c[2, 2]+1/4*c[3, 3],  

a[2] = 1/4*c[4, 3]+1/4*c[2, 1]-1/4*c[1, 2]-1/4*c[3, 4],  

a[3] = 1/4*c[2, 4]+1/4*c[3, 1]-1/4*c[1, 3]-1/4*c[4, 2],  

a[4] = -1/4*c[1, 4]-1/4*c[2, 3]+1/4*c[4, 1]+1/4*c[3, 2],  

a[5] = 1/4*c[4, 3]-1/4*c[2, 1]+1/4*c[1, 2]-1/4*c[3, 4],  

a[6] = 1/4*c[2, 4]-1/4*c[3, 1]+1/4*c[1, 3]-1/4*c[4, 2],  

a[7] = 1/4*c[1, 4]-1/4*c[2, 3]-1/4*c[4, 1]+1/4*c[3, 2],  

a[8] = 1/4*c[1, 1]-1/4*c[4, 4]+1/4*c[2, 2]-1/4*c[3, 3],  

a[9] = 1/4*c[3, 2]-1/4*c[1, 4]+1/4*c[2, 3]-1/4*c[4, 1],  

a[10] = 1/4*c[2, 4]+1/4*c[3, 1]+1/4*c[1, 3]+1/4*c[4, 2],  

a[11] = 1/4*c[1, 4]+1/4*c[2, 3]+1/4*c[4, 1]+1/4*c[3, 2],  

a[12] = 1/4*c[1, 1]-1/4*c[4, 4]-1/4*c[2, 2]+1/4*c[3, 3],  

a[13] = 1/4*c[4, 3]-1/4*c[2, 1]-1/4*c[1, 2]+1/4*c[3, 4],  

a[14] = 1/4*c[2, 4]-1/4*c[3, 1]-1/4*c[1, 3]+1/4*c[4, 2],  

a[15] = 1/4*c[4, 3]+1/4*c[2, 1]+1/4*c[1, 2]+1/4*c[3, 4],  

a[16] = 1/4*c[1, 1]+1/4*c[4, 4]-1/4*c[2, 2]-1/4*c[3, 3],  

>
 

>