The commutation relations for the quantum mechanical angular momentum operators are often expressed as L x L = i hbar L where the x is the vector cross product and hbar is Planck's constant divided by 2 pi and L is the operator. Setting up vector operators turned out to be complicated than I expected. I had to set up vectors of functions, rather than functions returning vectors, so I could access the components as functions. Suggestions for improvement of the code would be appreciated Note that I have a vector operator dot product function defined as well as the cross product, but that is not used in the proof. .

 Position operator 

In[1]:= xop = x*# &; yop = y*# &; zop = z*# &;

In[2]:= rop = {xop, yop, zop};

 Momentum operator

In[3]:= pxop = -I \[HBar] D[#, x] &; pyop = -I \[HBar] D[#, 
    y] &; pzop = -I \[HBar] D[#, z] &;

In[4]:= pop = {pxop, pyop, pzop};

 Operator Dot Product Function

In[5]:= opdot[{a_, b_, c_}, {d_, e_, 
   f_}] := (Composition[a, d][#] + Composition[b, e][#] + 
    Composition[c, f][#]) &

 Operator Cross Product Function

In[6]:= opcross[{a_, b_, c_}, {d_, e_, 
   f_}] := {(-Composition[c, e][#] + 
     Composition[b, f][#]) &, (Composition[c, d][#] - 
     Composition[a, f][#]) &, (-Composition[b, d][#] + 
     Composition[a, e][#]) &}

 Angular Momentum operator L

In[7]:= lop = opcross[rop, pop];

 L^2

In[8]:= lsqrop = opdot[lop, lop];

 L x L

In[9]:= lxlop = opcross[lop, lop];

 L x L acting on an arbitary function

In[10]:= lxlopf = Through[lxlop[f[x, y, z]]] // FullSimplify;

i \[HBar] L acting on an arbitrary function

In[11]:= lopf = I \[HBar]*Through[lop[f[x, y, z]]] // FullSimplify;

 L x L = i \[HBar] L

In[12]:= Thread[lxlopf == lopf]

Out[12]= True

Attachments: