How to speed up my calculation that involves integral2 ?

To start, here is a portion of my code that utilizes for loops:

N = 128;
xmax = 1.1*10^(-3);
ymax = 1.1*10^(-3);

for n1 = 1:N
    xprime = -1*xmax  + 2*(n1-1)*(xmax/(N-1));
    for n2 = 1:N
        yprime = -1*ymax + 2*(n2-1)*(ymax/(N-1));

          kernel = @(x,y)(wavenumGRIN/(1i.*2.*pi.*B))*exp(1i.*wavenumGRIN.*n0_GRIN.*L)*exp(((1i.*wavenumGRIN)/(2.*B)).*((xprime.^2 + yprime.^2).*D + (x.^2 + y.^2).*A - 2.*(x.*xprime + y.*yprime))).*exp(-1*(x.^12 + y.^12)/((200e-6)^12));
          wave_out = integral2(kernel, xminGRIN, xmaxGRIN, yminGRIN, ymaxGRIN, 'Method','iterated','AbsTol',0.001,'RelTol',0.001);
          GRINoutput(n1,n2) = wave_out;
          xprime2(n1,n2) = xprime;
          yprime2(n1,n2) = yprime;

      end
  end

scaledGRINout = GRINoutput/(abs(max(max(GRINoutput))));

I tried using the tiled method, which is faster than iterated, but unfortunately I received this error message: "Warning: Reached the maximum number of function evaluations (10000). The result passes the global error test."

I tried relaxing the tolerances, but that did not help much and I would prefer to maintain my current tolerances.

Would vectorizing my code work/help? How can I do this?

Thanks!