Actually, I want do numerical integration, I use the Simpson method, the integration in on k and w, so there are two loops. I use Table to replace the for loop, but the speed almost the same. I also write the Matlab code, it's much fast. this result is right, I want to speed up it , how can I do?

unispec[w_, n_, d_, min_, max_, epsi10_, epsi1e_, epsi20_, epsi2e_] :=
  Module[{},
  dkx = (max - min)/n;
  kx = Table[0, {i, 1, n + 1}];
  swkx = Table[0, {i, 1, n + 1}];
  k0 = w/c0;
  Table[{kx[[ind]] = min + (ind - 1)*dkx;
    kz10 = 
     If[Im[Sqrt[epsi10 k0^2 - kx[[ind]]^2]] >= 0, 
      Sqrt[epsi10 k0^2 - kx[[ind]]^2], -Sqrt[
        epsi10 k0^2 - kx[[ind]]^2]];
    kz20 = 
     If[Im[Sqrt[epsi20 k0^2 - kx[[ind]]^2]] >= 0, 
      Sqrt[epsi20 k0^2 - kx[[ind]]^2], -Sqrt[
        epsi20 k0^2 - kx[[ind]]^2]];
    kz1e = 
     If[Im[Sqrt[epsi10 k0^2 - kx[[ind]]^2 epsi10/epsi1e]] >= 0, 
      Sqrt[epsi10 k0^2 - kx[[ind]]^2 epsi10/epsi1e],
      -Sqrt[epsi10 k0^2 - kx[[ind]]^2 epsi10/epsi1e]];
    kz2e = 
     If[Im[Sqrt[epsi20 k0^2 - kx[[ind]]^2 epsi20/epsi2e]] >= 0, 
      Sqrt[epsi20 k0^2 - kx[[ind]]^2 epsi20/epsi2e],
      -Sqrt[epsi20 k0^2 - kx[[ind]]^2 epsi20/epsi2e]];
    kz3 = Sqrt[k0^2 - kx[[ind]]^2];
    rs31 = (kz3 - kz10)/(kz3 + kz10);
    rp31 = (epsi10 kz3 - kz1e)/(epsi10 kz3 + kz1e);
    rs32 = (kz3 - kz20)/(kz3 + kz20);
    rp32 = (epsi20 kz3 - kz2e)/(epsi20 kz3 + kz2e);
    tpstemp = Abs[1 - rs31 rs32 Exp[2 I kz3 d]];
    tpptemp = Abs[1 - rp31 rp32 Exp[2 I kz3 d]];
    etemp = Exp[-2 Im[kz3] d];
    If[kx[[ind]] <= 
      k0, {trans = (1 - Abs[rs31]^2) (1 - Abs[rs32]^2)/
          tpstemp^2 + (1 - Abs[rp31]^2) (1 - Abs[rp32]^2)/tpptemp^2;
      swkx[[ind]] = kx[[ind]] trans/(4 Pi^2);},
     {trans = 
       4 etemp Im[rs31] Im[rs32]/Abs[1 - rs31 rs32 etemp]^2 + 
        4 etemp Im[rp31] Im[rp32]/Abs[1 - rp31 rp32 etemp]^2;
      swkx[[ind]] = kx[[ind]] trans/(4 Pi^2)}];}, {ind, 1, n + 1}];
  swtot = (swkx[[1]] + 4 Sum[swkx[[k]], {k, 2, n, 2}] + 
      2 Sum[swkx[[k]], {k, 3, n - 1, 2}] + swkx[[n + 1]])*dkx/3]


Qtd = {};
For[i = 1, i <= 3, i = i + 0.2,
 d = 10^i nm;
 Print[10^i];
 Print[DateString[]];
 Table[
  {wi = wstart + (ii - 1) dw;
   a1 = (wi - 10)/c0;
   ae1 = (wi + 10)/c0;
   ae2 = Max[0.6/d, const1 wi/c0];
   ae3 = Max[5/d, const2 wi/c0];
   Q1[[ii]] = h wi/(Exp[h wi/(kb T1)] - 1);
   Q2[[ii]] = h wi/(Exp[h wi/(kb T2)] - 1);
   twp[[ii]] = 
    unispec[wi, nk, d, a0, a1, epsi01[wi], epsie1[wi], epsi02[wi], 
     epsie2[wi]];
   twe1[[ii]] = 
    unispec[wi, nk, d, ae1, ae2, epsi01[wi], epsie1[wi], epsi02[wi], 
     epsie2[wi]];
   twe2[[ii]] = 
    unispec[wi, nk, d, ae2, ae3, epsi01[wi], epsie1[wi], epsi02[wi], 
     epsie2[wi]];
   sw[[ii]] = (twp[[ii]] + twe1[[ii]] + twe2[[ii]]) (Q1[[ii]] - 
       Q2[[ii]])}
  , {ii, 1, nw}]; 
 Qt = (sw[[1]] + 4 Sum[sw[[k]], {k, 2, nw - 1, 2}] + 
     2 Sum[sw[[k]], {k, 3, nw - 1, 2}] + sw[[nw]]) dw/3;
 AppendTo[Qtd, Qt];
 Print[DateString[]];
 ]

Attachments: