Let's look at the number of iterations you perform in each loop. Pulling out appropriate sections of your code to define the variables used to create the data over which you iterate (leaving it commented out because we don't actually need to run it, just to refer to it):
%{
twidth = 500e-15; % Time window [s]
%twidth = 30e-12;
n = 2^14; % Number of points
lambda0 = 795e-9; % Center wavelength [m]
W0 = 2*pi*c/lambda0; % Center frequency [rad/s]
FWHM = 34e-15; % Initial pulse duration [s]
Epsilon_0 = 8.8541878128*1e-12;
% Time domain
T = linspace(-twidth/2, twidth/2, n); % Time grid
% Spatial domain
x = linspace(-5e-3, 5e-3, 256);
E_x_y_t = zeros(length(x), length(x), length(T));
for i = 1:length(x)
for j = 1:length(x)
for k = 1:length(T)
E_x_y_t(i, j, k) = (1/Normalization)*gaussian_space(i, j) * fin_t(k);
end
end
end
%}
So your outer loop runs 256 times as does your middle loop. Your innermost loop runs 2^14 times. That's:
format longg
numInnermostExecutions = 256*256*(2^14)
just over a billion executions of the innermost loop. You haven't showed us what your fin_t function is doing, but I'd execute that function once on the T vector (if it's vectorized, or once per element of T if it's not) and then compute E_x_y_t using implicit expansion. No loops required. As an example of a similar operation:
M = magic(4)
fin_t_vector = (1:3).^2 % Let's say your function was fin_t = @(t) t.^2
E = M .* reshape(fin_t_vector, 1, 1, [])
I haven't looked too closely at the later loops but I suspect you may be able to use a similar approach to avoid executing your functions billions of times.
See this documentation page for more information on implicit expansion and how it may allow you to avoid performing scalar operations over and over and over again.
