https://mathematica.stackexchange.com/questions/156900/can-any-one-help-me-
make-my-program-work-faster

avec = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
nvec = {1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6};
svec = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
ne = 3;
 nμ = 12;
 δ = -50;
 β = 1;
 kdfxn[i_, j_] := If[i == j, 1, 0]

 cfxn = Block[{n1a, n1b, n2a, n2b},
  With[{code = Which[
   n1a == n1b && n2a == n2b,
   Evaluate[N[1/6 (1 - 3/(n1a^2 π^2) - 3/(n2a^2 π^2))]],

   n1a == n1b && n2a != n2b,
   Evaluate[N[(4 (1 + (-1)^(n2a + n2b)) n2a n2b)/((n2a^2 - n2b^2)^2 π^2)]],

   n1a != n1b && n2a == n2b,
   Evaluate[N[(4 (1 + (-1)^(n1a + n1b)) n1a n1b)/((n1a^2 - n1b^2)^2 π^2)]],

   True,
   Evaluate[N[-((32 (-1 + (-1)^(n1a + n1b)) (-1 + (-1)^(n2a + n2b)) n1a n1b
   n2a n2b)/((n1a^2 - n1b^2)^2 (n2a^2 - n2b^2)^2 π^4))]]
   ]},
   Compile[{{n1a, _Integer}, {n1b, _Integer}, {n2a, _Integer}, {n2b,
   _Integer}},
    code,
  CompilationTarget -> "C"
   ]
    ]];

chmat = With[{ccfxn = cfxn, kkdfxn = kdfxn},
   Compile[{{nm, _Integer}, {ne, _Integer}, {b, _Real}, {d, _Real},
   {avec, _Real, 2}, {nvec, _Real, 1}, {svec, _Real, 1}},
   Block[{sz0, sz1, sz2, sz3, n0, n1, n2, n3, h1, h2, tmp, tmp2,
   tmp21, kf01, kf23, kf13, kf02},
   Table[n0 = Compile`GetElement[nvec, nm0];
   n1 = Compile`GetElement[nvec, nm1];
   sz0 = Compile`GetElement[svec, nm0];
   sz1 = Compile`GetElement[svec, nm1];
   tmp = 0.;
   Do[sz2 = Compile`GetElement[svec, nm2];
   sz3 = Compile`GetElement[svec, nm3];
   n2 = Compile`GetElement[nvec, nm2];
   n3 = Compile`GetElement[nvec, nm3];
   tmp2 = ccfxn[n1, n0, n3, n2];
   tmp21 = ccfxn[n1, n3, n0, n2];

   kf01 = kkdfxn[sz0, sz1];
   kf23 = kkdfxn[sz2, sz3];
   kf13 = kkdfxn[sz1, sz3];
   kf02 = kkdfxn[sz0, sz2];

   Do[
    tmp += (tmp2 kf23 kf01 - tmp21 kf13 kf02) Compile`GetElement[
       avec, j, nm3] Compile`GetElement[avec, j, nm2], {j, 1,
     ne}], {nm2, 1, nm}, {nm3, 1, nm}];
    d tmp +
    If[nm0 == nm1, (n0^2 Pi^2 + b Compile`GetElement[svec, nm0]),
    0.], {nm0, 1, nm}, {nm1, 1, nm}]], CompilationTarget -> "C",
  CompilationOptions -> {"InlineCompiledFunctions" -> True},
  RuntimeOptions -> "Speed"]];

Table[
  hmat = chmat[n[Mu], ne, [Delta], [Beta], avec, nvec, svec];
  {evals, evecs} = Eigensystem[hmat];
  pos = Ordering[evals][[1 ;; ne]];
  bvec = Map[x [Function] If[Total[x] < 0, -x, x], evecs[[pos]]];
  residual = Max[Abs[avec - bvec]];
  avec = bvec;
  {residual, Total[evals[[pos]]]},
  {j, 1, 30}]