When I define a Piecewise function containing a solution to a partial differential equation, the function cannot be evaluated numerically

I need to find an integral, in which there is a function

YO [x_] = 
  Sum[E^(((I h PI)/(λ 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 the domain of YO 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 get the 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 g[xo, zo] by FindArgMax, Whether it is the function of YO before expansion or the function of YO after integration, It will report an error that the function value is not a real number. But my function value is the square after taking the modulus, 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 code.

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, and it is proved that the solution is completely correct.

ClearAll["Global`*"];
θ = 0;
λ = 1.2398/(19.5 10^3);
f = 4.72 10^3;
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][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 π)/(λ f) x^2) Subscript[Y, h][x, z = 13.5],
        {h, -5, -1}],
      20 >= x >= 0}, 
     {Sum[
        E^((I h π)/(λ 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 of g[xo, zo] 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 λ 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 NIntegrate and FindArgMax will give me an error. The message from NIntegrate tells me that some values are not numerical, and the message from FindArgMax tells me that some values 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.

This is my unexpanded YO[x], which is defined over {x, 0, 30}. With this function, I can run YO[3] to calculate the specific value, but if I expand the domain to {x, - 30,30}, which is the piecewise function in the previous code, I can only get itself by running YO[3].

YO[x_]=
  Sum[E^(((I h Pi)x^2)/(λ f)) Subscript[Y, h][x, z = 13.5], 
  {h, -5, -1}] /. s; 

I wonder if this is the reason why plotting after integration produces complaints about non-numerical values, but no matter which YO[x] I use to calculate the maximum of g[xo,zo], it always complains about non-numerical values.

θ = 0;
λ = 1.2398/(19.5 10^3);
f = 4.72 10^3;
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, 30}, {z, 0, 30},
    Method ->
      {"MethodOfLines",
       "SpatialDiscretization" ->
         {"TensorProductGrid", "MaxPoints" -> 100,
          "MinPoints" -> 100, "DifferenceOrder" -> 4},
       "TemporalVariable" -> z}];

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

YO /@ {-20, -3, 3, 20} // Column