MatLab code speed and optimization. How to improve?

Question

Could someone please run this for me and tell me how long it takes for you? It took my laptop 60s. I can't tell if it's my laptop that's crappy or my code. Probably both. I just started learning MatLab, so I'm not yet familiar with which functions are better than others for specific tasks. If you have any suggestions on how I could improve this code, it would be greatly appreciated.

function gbp
clear; clc;

zi = 0;                                  % initial position
zf = 100;                                % final position

Ei = 1;                                  % initial electric field
c = 3*10^8;                              % speed of light
epsilon = 8.86*10^-12;                   % permittivity of free space
lambda = 1064*10^-9;                     % wavelength
k = 2*pi/lambda;                         % wave number
wi = 1.78*10^-3;                         % initial waist width (minimum spot size)
zr = (pi*wi^2)/lambda;                   % Rayleigh range
Ri = zi + zr^2/zi;                       % initial radius of curvature
qi = 1/(1/Ri-1i*lambda/(pi*wi^2));       % initial complex beam parameter
Psii = atan(real(qi)/imag(qi));          % Gouy phase
mat = [1 zf; 0 1];                       % transformation matrix

A = mat(1,1); B = mat(1,2); C = mat(2,1); D = mat(2,2);
qf = (A*qi + B)/(C*qi + D);              % final complex beam parameter
wf = sqrt(-lambda/pi*(1/imag(1/qf)));    % final spot size
Rf = 1/real(1/qf);                       % final radius of curvature
Psif = atan(real(qf)/imag(qf));          % final Gouy phase

% Hermite - Gaussian modes function
u = @(z, x, n, w, R, Psi) (2/pi)^(1/4)*sqrt(exp(1i*(2*n+1)*Psi)/(2^n*factorial(n)*w))*...
            hermiteH(n,sqrt(2)*x/w).*exp(-x.^2*(1/w^2+1i*k/(2*R))-1i*k*z);

% Complex amplitude coefficients function
a = @(n) exp(1i*k*zi)*integral(@(x) Ei.*conj(u(zi, x, n, wi, Ri, Psii)),-2*wi,2*wi);

%----------------------------------------------------------------------------
xlisti = -0.1:1/10000:0.1;              % initial x-axis range
xlistf = -0.1:1/10000:0.1;              % final x-axis range
nlist = 0:2:20;                         % modes range

     function Eiplot
        Efieldi = zeros(size(xlisti));
        for nr = nlist
            Efieldi = Efieldi + a(nr).*u(zi, xlisti, nr, wi, Ri, Psii)*exp(-1i*k*zi);
        end
        Ii = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldi);
    end

    function Efplot
        Efieldf = zeros(size(xlistf));
        for nr = nlist
            Efieldf = Efieldf + a(nr).*u(zf, xlistf, nr, wf, Rf, Psif)*exp(-1i*k*zf);
        end
        If = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldf);
    end

Eiplot
Efplot

plot(xlisti,real(Ii),xlistf,real(If))

xlabel('x(m)')                       % x-axis label
ylabel('I(W/m^2)')                   % y-axis label
end

Answer 1

The cost is coming from the calls to hermiteH -- for every call, this creates a new function using symbolic variables, then evaluates the function at your input. The key to speeding this up is to pre-compute the hermite polynomial functions then evaluate those rather than create them from scratch each time (speedup from ~26 seconds to around 0.75 secs on my computer).

With the changes:

function gbp

x = sym('x');

zi = 0;                                  % initial position
zf = 100;                                % final position

Ei = 1;                                  % initial electric field
c = 3*10^8;                              % speed of light
epsilon = 8.86*10^-12;                   % permittivity of free space
lambda = 1064*10^-9;                     % wavelength
k = 2*pi/lambda;                         % wave number
wi = 1.78*10^-3;                         % initial waist width (minimum spot size)
zr = (pi*wi^2)/lambda;                   % Rayleigh range
Ri = zi + zr^2/zi;                       % initial radius of curvature
qi = 1/(1/Ri-1i*lambda/(pi*wi^2));       % initial complex beam parameter
Psii = atan(real(qi)/imag(qi));          % Gouy phase
mat = [1 zf; 0 1];                       % transformation matrix

A = mat(1,1); B = mat(1,2); C = mat(2,1); D = mat(2,2);
qf = (A*qi + B)/(C*qi + D);              % final complex beam parameter
wf = sqrt(-lambda/pi*(1/imag(1/qf)));    % final spot size
Rf = 1/real(1/qf);                       % final radius of curvature
Psif = atan(real(qf)/imag(qf));          % final Gouy phase


% Hermite - Gaussian modes function
nlist = 0:2:20;     % modes range

% precompute hermite polynomials for nlist
hermites = {};
for n = nlist
    if n == 0
        hermites{n + 1} = @(x)1.0;
    else
        hermites{n + 1} = matlabFunction(hermiteH(n, x));
    end
end

u = @(z, x, n, w, R, Psi) (2/pi)^(1/4)*sqrt(exp(1i*(2*n+1)*Psi)/(2^n*factorial(n)*w))*...
            hermites{n + 1}(sqrt(2)*x/w).*exp(-x.^2*(1/w^2+1i*k/(2*R))-1i*k*z);

% Complex amplitude coefficients function
a = @(n) exp(1i*k*zi)*integral(@(x) Ei.*conj(u(zi, x, n, wi, Ri, Psii)),-2*wi,2*wi);

%----------------------------------------------------------------------------
xlisti = -0.1:1/10000:0.1;              % initial x-axis range
xlistf = -0.1:1/10000:0.1;              % final x-axis range

     function Eiplot
        Efieldi = zeros(size(xlisti));
        for nr = nlist
            Efieldi = Efieldi + a(nr).*u(zi, xlisti, nr, wi, Ri, Psii)*exp(-1i*k*zi);
        end
        Ii = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldi);
    end


    function Efplot
        Efieldf = zeros(size(xlistf));
        for nr = nlist
            Efieldf = Efieldf + a(nr).*u(zf, xlistf, nr, wf, Rf, Psif)*exp(-1i*k*zf);
        end
        If = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldf);
    end

Eiplot
Efplot

plot(xlisti,real(Ii),xlistf,real(If))

xlabel('x(m)')                       % x-axis label
ylabel('I(W/m^2)')                   % y-axis label
end