I need to find an integral, in which there is a function YO(x), YO (x_)= sum (e ^ (((I * h * PI) / ( (lambda) * f)) * x ^ 2) * subscript (Y, h) (x, z = 13.5), {h, – 5, – 1}) /. S (where subscript (Y, h) (x, z) is a solution of partial differential equation). I run YO (3) in the positive interval of {x, 0, 30} and find that I can get a numerical value. After that, I want to extend YO which only contains x positive direction {x, 0, 30} to {x, -30,30}. The value of negative interval upper function and positive interval are symmetric. So I manually defined a piecewise function, and tried it with Piecewise and If functions. When the function YO (x) is in the negative interval, the value is equal to YO (Abs (x)), but I found that after calculation, running YO (3) can not get the numerical value, so I can only return YO (3) itself.

Based on the above situation, it is OK for me to integrate with the non expanded YO(x), but if I integrate with the expanded YO(x) with Piecewise or If function, I will report an error: the sampling point is non numerical at a certain point. At the same time, It also report errors when I find the maximum value of the integral result by Findargmax: Whether it is the function of YO before expansion or the function of YO after integration, It will report an error The display function value is not a real number. But my function value is the square after taking the module, which must be a real number. The relevant knowledge that I have been looking for for a long time these days has not been solved, so I would like to know what’s wrong with my program.

In addition, I can get the specific value by directly integrating the function YO (x) which is not expanded, but I can’t get the numerical value by integrating the function YO (x) which is expanded by piecewise and if, so I can only get the function formula which contains subscript (Y, h)(x,z).

This is my program for solving partial differential equations

```
ClearAll("Global`*");
(Theta) = 0;
(Lambda) = 1.2398/(19.5 10^3);
f = 4.72 10^3;
Subscript((Delta), 1) = 1.274/10^6;
Subscript((Delta), 2) = 4.304/10^6;
Subscript(bt, 1) = 5.254/10^9;
Subscript(bt, 2) = 2.435/10^7;
(Chi)1 = -2 Subscript((Delta), 1) + 2 I Subscript(bt, 1);
(Chi)2 = -2 Subscript((Delta), 2) + 2 I Subscript(bt, 2);
(CapitalDelta)(Chi) = (Chi)1 - (Chi)2;
k = 2 ((Pi)/(Lambda));
Table((Beta)(b) = -((2 (b x) Sin((Theta)))/f) - ((b x)/f)^2, {b, -5,
5});
Table((Chi)(
a) = ((CapitalDelta)(Chi) (1 - (-1)^
Abs(a)))/(2 I a (Pi)), {a, -10, -1});
Table((Chi)(
a) = ((CapitalDelta)(Chi) (1 - (-1)^Abs(a)))/(2 I a (Pi)), {a,
1, 10});
Table((Chi)(a) = ((Chi)1 + (Chi)2)/2, {a, 0, 0});
eqns = Join(Table((Sin((Theta)) + (h x)/f) !(
*SubscriptBox(((PartialD)), (x))(
(*SubscriptBox((Y), (h)))(x, z))) + Cos((Theta)) !(
*SubscriptBox(((PartialD)), (z))(
(*SubscriptBox((Y), (h)))(x,
z))) == ((I (Pi)) ((Beta)(h) Subscript(Y, h)(x, z) + !(
*UnderoverscriptBox(((Sum)), (l = (-5)), (5))((Chi)(
h - l)
(*SubscriptBox((Y), (l)))(x, z)))))/(Lambda), {h, -5, 5}),
Table(Subscript(Y, h)(x, 0) == If(h == 0, 1, 0), {h, -5, 5}));
s = NDSolve(eqns,
Table(Subscript(Y, h)(x, z), {h, -5, 5}), {x, 0, 30}, {z, 0, 30},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 100, "MinPoints" -> 100,
"DifferenceOrder" -> 4}, "TemporalVariable" -> z});
```

This is my program to define YO(x)

```
YO(x_?NumericQ) :=
Piecewise({{Sum(
E^((I h (Pi))/((Lambda) f) x^2)
Subscript(Y, h)(x, z = 13.5), {h, -5, -1}),
20 >= x >= 0}, {Sum(
E^((I h (Pi))/((Lambda) f) x^2)
Subscript(Y, h)((Abs(x)), z = 13.5), {h, -5, -1}),
0 > x >= -20}}, {x, -20, 20}) /. s
```

And this is what I want to do: integrate, find the maximum and plot the density around zo at the maximum

```
f1(xo_?NumericQ, zo_?NumericQ) :=
NIntegrate((E^(I k zo) YO(x) E^(((I k) x^2)/(2 zo))
E^(((I k) xo^2)/(2 zo))
E^(-(((I k) x xo)/zo)))/(I (Lambda) zo^0.5), {x, -20,
20}) /. s;
g(xo_?NumericQ, zo_?NumericQ) := Abs(f1(xo, zo)^2);
FindArgMax({g(xo, zo), -20 <= xo <= 20 && 3000 <= zo <= 5000}, {{xo,
0}, {zo, 4000}})
DensityPlot(g(xo, zo), {xo, -20, 20}, {zo, 3000, 5000},
PlotLegends -> Automatic)
```

But in fact, both the integral and the maximum value will prompt me to report an error. The integral will prompt me that some places are not numerical, and the maximum value will prompt me that some places are not real. But it’s obvious that the absolute value followed by the square must be real and numerical.

In addition, when I integrate with the unexpanded YO(x), I can integrate and draw a graph, but I never succeed in finding the maximum value.