I am constructing a matrix of the form
$$
T_{i_aj_ak_al_a,i_bj_bk_bl_b} = T_{i_ai_b}\delta_{j_aj_b}\delta_{k_ak_b}\delta_{l_al_b} + T_{j_aj_b}\delta_{i_ai_b}\delta_{k_ak_b}\delta_{l_al_b} + T_{k_ak_b}\delta_{i_ai_b}\delta_{j_aj_b}\delta_{l_al_b} + T_{l_al_b}\delta_{i_ai_b}\delta_{j_aj_b}\delta_{k_ak_b}
$$

with

$$
T_{i_ai_b} = \frac{\hbar (-1)^{i_a-i_b}}{2m_w \Delta w^2} \left\{ \begin{array}{ll} \dfrac{\pi^2}{3},  & i_a = i_b \\ \dfrac{2}{(i_a-i_b)^2}, & i_a \neq i_b \end{array} \right\}
$$

I did this in *Mathematica* with the following code:

    Twxyz =
      SparseArray[
       Flatten[
        Table[
         If[ia == ib,
            \[Pi]^2/3,
            2/(ia - ib)^2]*(hbar^2*(-1)^(ia - ib))/(2 mw *dw^2)*
           KroneckerDelta[ja, jb]*KroneckerDelta[ka, kb]*
           KroneckerDelta[la, lb]
          +
          If[ja == jb,
            \[Pi]^2/3,
            2/(ja - jb)^2]*(hbar^2*(-1)^(ja - jb))/(2 mx*dx^2)*
           KroneckerDelta[ia, ib]*KroneckerDelta[ka, kb]*
           KroneckerDelta[la, lb]
          +
          If[ka == kb,
            \[Pi]^2/3,
            2/(ka - kb)^2]*(hbar^2*(-1)^(ka - kb))/(2 my*dy^2)*
           KroneckerDelta[ia, ib]*KroneckerDelta[ja, jb]*
           KroneckerDelta[la, lb]
          +
          If[la == lb,
            \[Pi]^2/3,
            2/(la - lb)^2]*(hbar^2*(-1)^(la - lb))/(2 mz*dz^2)*
           KroneckerDelta[ia, ib]*KroneckerDelta[ja, jb]*
           KroneckerDelta[ka, kb]
         ,
         {lb, npointsgrid},
         {la, npointsgrid},
         {kb, npointsgrid},
         {ka, npointsgrid},
         {jb, npointsgrid},
         {ja, npointsgrid},
         {ib, npointsgrid},
         {ia, npointsgrid}]
        , {{1, 3, 5, 7}, {2, 4, 6, 8}}]
       ];

This approach works ok. The problem is that it takes several hours to build the matrix due to the many if statements it must process for each matrix element. Any suggestions on how to speed up the matrix construction?

Many thanks in advance!

EDIT:

The table is numeric.

Typical values for variables are

    hbar = 6.62607004*^-34/(2*Pi)
    npointsgrid = 13
    mx*dx^2 = 3.76547*^-49

In the end, the matrix is constructed such that $$ m_w\Delta w^2 = m_x\Delta x^2 = m_y \Delta y^2 = m_z \Delta z^2$$  

    In[307]:= mw*dw^2
              mx*dx^2
              my*dy^2
              mz*dz^2

    Out[307]= 3.76547*10^-49

    Out[308]= 3.76547*10^-49

    Out[309]= 3.76547*10^-49

    Out[310]= 3.76547*10^-49