My current solution involving a very disgusting for loop:

clf
clear

# initiation
syms x y z lambda pi;
SUMMATION = 0;

# numeric meshgrid specification
N = 5;
start_value = ((N-1)/2) * (sym(-136)/100000000)
end_value = ((N-1)/2) * (sym(136)/100000000)

# numeric meshgrid generation
xi = linspace(start_value,end_value,N);
eta = linspace(start_value,end_value,N);

[XI_numeric,ETA_numeric] = meshgrid(xi,eta)

# symbolic meshgrid generation
XI_symbolic = sym("xi",[N N]);
ETA_symbolic = sym("eta",[N N]);

[XI_symbolic_rows, XI_symbolic_cols] = size(XI_symbolic);

# iterative summation
for I = 1:XI_symbolic_rows
  for J = 1:XI_symbolic_cols

    element_symbolic = exp(-2*pi*1i * ( (x*XI_symbolic(I,J))/(lambda*z) + (y*ETA_symbolic(I,J))/(lambda*z) ));
    element_numeric = subs(element_symbolic , {XI_symbolic(I,J),ETA_symbolic(I,J)} , {XI_numeric(I,J),ETA_numeric(I,J)})

    SUMMATION = SUMMATION + element_numeric;

  end
end

disp(SUMMATION)