# A question about Plot an NIntegrate and FindMaximum is too slow

I need to DensityPlot an NIntegrate where near the maximum piont. So I need to use FindMaximumto find the my plot range, but it has been running for more than 48 hours and has not given any results if I don't define starting points and intervals, if I define the starting point and the interval, then the time will decrease as the interval decreases. So I think my integral function needs to speed up, I used Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0} method, it does improve the speed, but when I run to find the maximum, it obviously gives a wrong answer and some errors, so I use "LocalAdaptive", But it made the program run for a few hours too.

g[xo,zo] is the function that I want the maximum value and plot the density diagram, Its physical meaning is the wave field intensity of light. It's defined by the following functions.

f1[xo_?NumericQ, zo_?NumericQ] :=
NIntegrate[(E^(I k zo)  E^(((I k) xo^2)/(2 zo))
YO[x] E^(((I k) x^2)/(2 zo)) E^(-(((I k) x xo)/zo)))/
Sqrt[I λ zo], {x, -30000, 30000},
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
MinRecursion -> 3, WorkingPrecision -> 10];

g[xo_?NumericQ, zo_?NumericQ] := Abs[f1[xo, zo]^2];

FindMaximum[{g[xo, zo], -300. <= xo <= 300. && 4000000. <= zo <= 5000000}, {{xo, 0}, {zo, 4720000}}]



The above program in the case that I know the maximum interval with Findmaximum can calculate the result in the time I can receive. It can output

{1634.28, {xo -> 2.19213*10^-10, zo -> 4.71831*10^6}}


It spend 3563.0903136 seconds. And DensityPlot takes Longer Longer time.

DensityPlot[g[xo, zo], {zo, 4698310, 4738310}, {xo, -200, 200},
PlotLegends -> Automatic, ColorFunction -> "Rainbow", Mesh -> None,
MaxRecursion -> 5, PlotPoints -> 30, PlotRange -> All]


But if I don't know the FindMaximum interval, I need to:

FindMaximum[{g[xo, zo], -300. <= xo <= 300. && 0. <= zo <= 5000000}, {{xo, 0}, {zo, 4720000}}]


And it took me dozens of hours and it kept running so I stopped the calculation.

So First:
I want to know if there is any NIntegrate method that allows me to find the maximum value accurately and within the time not too long, even if it's a little bit less accurate.

Second:
If there are some ways to make DensityPlot faster? I tried some other ways to integrate and Plot quickly. It took me some times, but I failed because I was not very good at mathematics.Like Henrik Schumacher's answer

The designed focus is zo= 4.72 mm = 4.72*10^6 nm, so the focus should be around (xo=0,zo=4720000) which makes me know the interval of FindMaximum .

Ps: what I did is the content of a document and I want to repeat it.the document link

The code definition of YO[x] is as follows, it contains a solution of the partial differential equation， and it is proved that the solution is completely correct.

Clear["Global*"]
t = AbsoluteTime[];
θ = 1.1/1000;
λ = 1.2398/19.5 ;
f = 4.72 10^6;
Subscript[δ, 1] = 1.274/10^6;
Subscript[δ, 2] = 4.304/10^6;
Subscript[bt, 1] = 5.254/10^9;
Subscript[bt, 2] = 2.435/10^7;
χ1 = -2 Subscript[δ, 1] + 2 I Subscript[bt, 1];
χ2 = -2 Subscript[δ, 2] + 2 I Subscript[bt, 2];
Δχ = χ1 - χ2;
k = 2 (π/λ);
Table[β[b] = -((2 (b x) Sin[θ])/f) - ((b x)/f)^2, {b, -5,5}];
Table[χ[a] = (Δχ (1-(-1)^Abs[a]))/(2 I a π), {a, -10, -1}];
Table[χ[a] = (Δχ (1-(-1)^Abs[a]))/(2 I a π), {a, 1, 10}];
Table[χ[a] = (χ1 + χ2)/2, {a, 0, 0}];
eqns =
Join[
Table[
(Sin[θ] + (h x)/f) D[Subscript[Y, h][x, z], x] +
Cos[θ] D[Subscript[Y, h][x, z], z] ==
((I Pi)
(β[h] Subscript[Y, h][x, z] +
Sum[χ[h - l] Subscript[Y, l][x, z], {l, -5, 5}]))/λ,{h, -5, 5}],
Table[Subscript[Y, h][x, 0] == If[h == 0, 1, 0], {h, -5, 5}]];
s = NDSolve[eqns,
Table[Subscript[Y, h], {h, -5, 5}], {x, 0, 30000}, {z, 0, 30000},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 100, "MinPoints" -> 100,
"DifferenceOrder" -> 4}, "TemporalVariable" -> z}];

Clear[YO]
YO[x_?NumericQ /; (0 > x >= -30000)] := YO[-x]
YO[x_?NumericQ] =
First[Sum[E^((I h π)/(λ f) x^2) Subscript[Y, h][x, 13500], {h, -5, -1}] /. s];


This is the picture should:

This is what the picture I have make, and I just calculate h=-1 of YO[x] for calculating {h,-5,-1} is too long to calculate.

• To @shrocat Haven't you seen the error message of NDSolve for defining the YO: NDSolve::bcart: Warning: An insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution. >>  . Therefore it makes no sense to calculate further. Give additional boundary conditions. Commented Jan 15, 2021 at 6:17
• @Akku14 It sure have the error message, but I think it is fine because using YO[x] I can calculate the same result as the literature, The problem is that's too slow, If I do difference integration useing other languages like Matlab to Plot, it reduced accuracy but increased speed. So I wonder if I can adjust the parameter of the Method of NIntegrate` to reduce the accuracy within allowable limits for speeding up calculation. Commented Jan 15, 2021 at 10:27