Use MATLAB's "int" twice:
syms a x1 y1 x y z lambda
f = exp(-(x1^2+x2^2)/a^2)*exp(1I*pi/(lambda*z)*((x-x1)^2+(y-y1)^2));
E = int(int(f,x1,-3*a,3*a),x2,-3*a,3*a)
Alternatively:
syms a x1 y1 x y z lambda
f1 = exp(-x1^2/a^2+1I*pi/(lambda*z)*(x-x1)^2);
f2 = exp(-y1^2/a^2+1I*pi/(lambda*z)*(y-y1)^2);
E = int(f1,x1,-3*a,3*a)*int(f2,x2,-3*a,3*a)
Best wishes
Torsten.
