In[]:=
$Version
Out[]=
14.1.0 for Microsoft Windows (64-bit) (July 16, 2024)
In[]:=
ClearAll["`*"];Remove["`*"];
In[]:=
testSquareF[x_?NumericQ,y_?NumericQ,w_?NumericQ,z_?NumericQ,λ_?NumericQ]:=NIntegrate
π8
4
z
8
4
z
+4
2
z
2
(y-η)
-
4
(y-η)
+48
6
z
-4
4
z
2
(y-η)
+3
2
z
4
(y-η)

2
(x-ξ)
+-8
4
z
+12
2
z
2
(y-η)
-15
4
(y-η)

4
(x-ξ)

32
7
z
λ

z
2
z
+
2
(-y+η)
+
2
(-x+ξ)
,{ξ,-w/2,w/2},{η,-w/2,w/2},PrecisionGoal->2;
With[{n=4(*Ifthenparameterissmaller,thanisbetterresolution*)},ListPlot3D[Partition[Flatten[ParallelTable[{x,y,Re[testSquareF[x,y,200,1000,1]]},{x,-150,150,n},{y,-150,150,n},Method->"FinestGrained"]],3],ColorFunction"Topographic",PlotLabelsHoldForm[Re[testSquareF[x,y]]]]]//AbsoluteTiming
Out[]=
578.993,

In[]:=
With[{n=4(*Ifthenparameterissmaller,thanisbetterresolution*)},ListPlot3D[Partition[Flatten[ParallelTable[{x,y,Im[testSquareF[x,y,200,1000,1]]},{x,-150,150,n},{y,-150,150,n},Method->"FinestGrained"]],3],ColorFunction"Topographic",PlotLabelsHoldForm[Im[testSquareF[x,y]]]]]//AbsoluteTiming
Out[]=
563.558,

(*---------------------------Betterspeedup------------------------------------------------------------------------------------------------*)
Code borrowed from here:
https://mathematica.stackexchange.com/questions/173253/mathematica-taking-forever-to-compute/173274#173274
ClearAll["`*"];Remove["`*"];
testSquareG=Block{x,y,w,z,λ,ξ,η},Withcode=N
π8
4
z
8
4
z
+4
2
z
2
(y-η)
-
4
(y-η)
+48
6
z
-4
4
z
2
(y-η)
+3
2
z
4
(y-η)

2
(x-ξ)
+-8
4
z
+12
2
z
2
(y-η)
-15
4
(y-η)

4
(x-ξ)

32
7
z
λ

z
2
z
+
2
(-y+η)
+
2
(-x+ξ)
,Compile[{{x,_Real},{y,_Real},{w,_Real},{z,_Real},{λ,_Real},{ξ,_Real},{η,_Real}},code,RuntimeAttributes{Listable},ParallelizationTrue];(*Canbemorespeedupadding:CompilationTarget->"C"ifyouhaveinstalledC-language.*)​​​​{pts,weights,errweights}=NIntegrate`GaussRuleData[9,$MachinePrecision];​​m=20;(*Ifthemparameterisgreater,thanisbetterresolution*)​​n=20;(*Ifthenparameterisgreater,thanisbetterresolution*)​​w=200.;(*w=200*)​​z=1000.;(*z=1000*)​​λ=1.;(*λ=1*)​​ξdata=Partition[Subdivide[-w/2,w/2,m],2,1];​​ηdata=Partition[Subdivide[-w/2,w/2,n],2,1];​​{ξ,η}=Transpose[Tuples[{Flatten[ξdata.{1.-pts,pts}],Flatten[ηdata.{1.-pts,pts}]}]];​​ω=Flatten[KroneckerProduct[Flatten[KroneckerProduct[Differences/@ξdata,weights]],Flatten[KroneckerProduct[Differences/@ηdata,weights]]]];W0={x,y}{x,y,Re[testSquareG[x,y,w,z,λ,ξ,η].ω]};(*RealPartoftestSquareGfunction!!!*)​​X=150.;(*Rangefor:X∈(-150,150)*)​​Y=150.;(*Rangefor:Y∈(-150,150)*)​​data0=Outer[W0,Subdivide[-X,X,100],Subdivide[-Y,Y,100]];//AbsoluteTiming//First​​W1={x,y}{x,y,Im[testSquareG[x,y,w,z,λ,ξ,η].ω]};(*ImaginaryPartoftestSquareGfunction!!!*)​​data1=Outer[W1,Subdivide[-X,X,100],Subdivide[-Y,Y,100]];//AbsoluteTiming//First
Out[]=
98.2391
Out[]=
106.061
In[]:=
ListPlot3D[Flatten[data0,1],PlotRangeAll,AxesLabel{"x","y","Z"},ColorFunction"Topographic",PlotLabelsHoldForm[Re[testSquareG[x,y]]]]​​ListPlot3D[Flatten[data1,1],PlotRangeAll,AxesLabel{"x","y","Z"},ColorFunction"Topographic",PlotLabelsHoldForm[Im[testSquareG[x,y]]]]
Out[]=
Out[]=