I managed to double the speed by simply moving some multiplications and divisions... AbsoluteTiming[ or2 = -2./Omegar^2.; oz2 = 2./Omegaz^2.; p = p^2.; Io*Table[ Exp[(p[[i, t, 1]] + p[[i, t, 2]])*or2 - p[[i, t, 3]]*oz2], {t, 1, Deltat + 1}, {i, n}]; ] A dot product is 30 to 40% faster still. AbsoluteTiming[ or2 = -2./Omegar^2.; oz2 = 2./Omegaz^2.; v = {or2, or2, -oz2}; p = p^2.; Io*Table[Exp[p[[i, t]].v],{t,1,Deltat+1},{i,n}]; ]