Strange result for an integral containing Meijer G-function

The question is about Mathematica's result for the integral
$$\int_0^\infty \frac{16}{\pi^2} \left( \frac{\sqrt{\pi}}{4} r^2 % G_{1,3}^{2,1}\left(\frac{r^2}{4}\Bigg| \begin{array}{l} 0 \\ 0,0,-\frac{1}{2} \\ \end{array} \right) % \right)^2 dr \,,$$
whereby the integrand will be denoted by $$f(r)$$ in the following. $$G$$ denotes the Meijer G-function (name of the corresponding Mathematica function: MeijerG). According to Mathematica's Integrate function this integral is $$1$$. However, with NIntegrate one obtains the following results:
$$\int_0^5 f(r) d r > 7 \,,$$ $$\int_0^{10} f(r) d r > 16 \,.$$ Furthermore Mathematica confirms that $$f(r) \geq 0 \,,\ \ \forall r \geq 0$$ holds. Using this together with the results from the numeric integrals one obtains $$\int_0^\infty f(r) d r = \int_0^{10} f(r) d r + \int_{10}^\infty f(r) d r > 16 + \int_{10}^\infty f(r) d r > 16 \,,$$ since the second integral has to be $$\geq 0$$ due to the non-negativity of the integrand. This contradicts the result from Mathematica's Integrate, which says that the integral is $$1$$.

The question is: Where is the fault leading to the contradiction? Which statement is true?
(I would naively assume that Integrate is wrong in this case.)

Note: In the case one doubts the results for the numeric integrals from NIntegrate, one can plot the integrand. These plots make the results seem very plausible. A screenshot of the plots and of the used Mathematica commands can be found below. Additionally, the used code is listed below for the case that some wants to test it.

Used Mathematica code

intgd[r_] := 16/\[Pi]^2*(1/4 Sqrt[\[Pi]] r^2 MeijerG[{{0}, {}}, {{0, 0}, {-(1/2)}}, r^2/4])^2
Integrate[intgd[r], {r, 0, Infinity}]
NIntegrate[intgd[r], {r, 0, 5}]
FullSimplify[intgd[r] >= 0, Assumptions -> {r >= 0}]
GraphicsRow[{Plot[intgd[r], {r, 0, 20}], Plot[intgd[r], {r, 0, 50}]}]

The used Mathematica version is 12.2.

Clear["Global`*"]

\$Version

(* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)

intgd[r_] :=
16/π^2*(1/4 Sqrt[π] r^2 MeijerG[{{0}, {}}, {{0, 0}, {-(1/2)}},
r^2/4])^2

lim = Limit[intgd[r], r -> Infinity]

(* 16/π^2 *)

For larger values of r, arbitrary-precision rather than machine precision is required for the Plot

Plot[intgd[r], {r, 0, 50},
WorkingPrecision -> 30,
PlotRange -> All,
GridLines -> {None, {lim}},
GridLinesStyle -> Directive[Red, AbsoluteThickness]] An alternate representation of the integrand for r > 0 is

intgd2[r_] = intgd[r] // FunctionExpand // Simplify[#, r > 0] &

(* (16 r^2 (CoshIntegral[r] Sinh[r] - Cosh[r] SinhIntegral[r])^2)/π^2 *)

intgd[r] == intgd2[r] // FullSimplify[#, r > 0] &

(* True *)

For this alternate representation, even greater precision is required in the Plot

Plot[intgd2[r], {r, 0, 50},
WorkingPrecision -> 50,
PlotRange -> All,
GridLines -> {None, {lim}},
GridLinesStyle -> Directive[Red, AbsoluteThickness]] As expected, the integral on {0, Infinity} does not converge

Integrate[intgd2[r], {r, 0, Infinity}] • So, this is a nice way to show that Integrate[intgd[r], {r, 0, Infinity}] produces an incorrect result. This implies that Mathematica's Integrate has a bug here, right? (And thanks for showing how to plot the integrand adequately!) Nov 5 '21 at 20:53