I'm trying to determine a general equation if possible. However, if we were to look at the actual value ranges, w will be between 100 and 250, x and y will be between -5w and +5w, z will be between w and 5000, and lambda will be between 0.5 and 1.5.

To make a 3D plot, first I would assign values to w, z, and lambda, something like

ts[x_,y_]:=testSquare/.{w->200,z->1000,lambda->1}

then plot as

Plot3D[ts[x,y],{x,-150,150},{y,-150,150}]

The integration over xi and nu would result in an equation in just x, y, z, w, and lambda. I did this using the second-order approximation and got accurate results for conditions that met the Fresnel diffraction requirements; here I'm looking for something where the second-order approximation results in errors exceeding 10%.

Change in the code to better resolution and better grid resolution:

   R = 2;(* If R parameter greater -grid resolution better,but the calculations take longer*)
   data0 = Outer[W0, Subdivide[-X, X, R*X], Subdivide[-Y, Y, R*Y]]; // AbsoluteTiming // First

The code is not mine, so for more information about it, please contact the creator Henrik Schumacher

Wow, I'm learning so much about speeding up plotting here...This version took 16.8 seconds to run on my computer, the previous best was nearly 150 seconds and my original was almost 400 seconds. At 16.8 seconds it isn't worth it any more to try to find a closed-form solution.

Thanks!

Thanks for the ideas. I am trying Assumptions; in the worst case I can try with numerical values for some of the variables. I did once take an integral that was solved after over 70,000 seconds (per AbsoluteTiming). It was also a little strange that the first time I ran this integration it gave me an error after 1716 seconds -- and the error was that I hadn't defined insideParams4.

I have found that making a table of values using ParallelTable, where I evaluated NIntegrate at 10,000 points, and it took significantly longer than using Plot3D. I haven't tried using Parallelize[NIntegrate]] because I have had many problems with that combination in the past.

The original integral has, for example, sqrt((xi-x)^2+(eta-y)^2+z^2). If I use a second-order approximation I get a closed-form solution but the integrand has errors of 10%. I might be able to accept 1%, which is why I'm trying this fourth-order approximation. I'm trying for a closed-form solution because calculating one for second order took about 1% as long as Plot3D on NIntegrate, and I was able to use Plot3D on a closed-form equation. But, as I said, second-order is not sufficient for this problem.

One thing I could also try is to compare the numerical integration time of the fourth-order approximation to that of the full equation using square roots. My instinct is that the square root would be equally fast or even faster, since the fourth-order version has three times as many variables.

Since I'm not a mathematician, I don't know much more than run this in Mathematica and see what happens. There are two standard approaches to this specific type of problem. One is the second-order approximation, which is insufficient for what I was trying to do; the other is limited to a size about 1/1000 of the area I need.