Problem plotting expression involving Generalized hypergeometric functions $_2F_2 left(.,.,. right)$

I’m trying to plot a graph for the following expectation

$$mathbb{E}left[ a mathcal{Q} left( sqrt{b } gamma right) right]=a 2^{-frac{kappa }{2}-1} b^{-frac{kappa }{2}} theta ^{-kappa } left(frac{, _2F_2left(frac{kappa }{2}+frac{1}{2},frac{kappa }{2};frac{1}{2},frac{kappa }{2}+1;frac{1}{2 b theta ^2}right)}{Gamma left(frac{kappa }{2}+1right)}-frac{kappa , _2F_2left(frac{kappa }{2}+frac{1}{2},frac{kappa }{2}+1;frac{3}{2},frac{kappa }{2}+frac{3}{2};frac{1}{2 b theta ^2}right)}{sqrt{2} sqrt{b} theta Gamma left(frac{kappa +3}{2}right)}right)$$
where $a$ and $b$ are constant values, $mathcal{Q}$ is the Gaussian Q-function, which is defined as $mathcal{Q}(x) = frac{1}{sqrt{2 pi}}int_{x}^{infty} e^{-u^2/2}du$ and $gamma$ is a random variable with Gamma distribition, i.e., $f_{gamma}(y) sim frac{1}{Gamma(kappa)theta^{kappa}} y^{kappa-1} e^{-y/theta} $ with $kappa > 0$ and $theta > 0$.

This equation was also found with Mathematica, so it seems to be correct. I’ve got the same plotting issue with Matlab.

Follows some examples, where I have checked the analytical results against the simulated ones.

When $kappa = 12.85$, $theta = 0.533397$, $a=3$ and $b = 1/5$ it returns the correct value $0.0218116$.

When $kappa = 12.85$, $theta = 0.475391$, $a=3$ and $b = 1/5$ it returns the correct value $0.0408816$.

When $kappa = 12.85$, $theta = 0.423692$, $a=3$ and $b = 1/5$ it returns the value $-1.49831$, which is negative. However, the correct result should be a value around $0.0585$.

When $kappa = 12.85$, $theta = 0.336551$, $a=3$ and $b = 1/5$ it returns the value $630902$. However, the correct result should be a value around $0.1277$.

Therefore, the issue happens as $theta$ decreases. For values of $theta > 0.423692$ the analytical matches the simulated results. The issue only happens when $theta <= 0.423692$.

I’d like to know if that is an accuracy issue or if I’m missing something here and if there is a way to correctly plot a graph that matches the simulation. Perhaps there is another way to derive the above equation with other functions or there might be a way to simplify it and get more accurate results.

prec = 100;
m = 16;
a = 4 (1 - (1/Sqrt[m]));
b = 3/(m - 1);

pdb = -12;
p = 10^(pdb/10);

betag = (1.5)^2;

betah = (2.3)^2;

bits = 8;

Q = 2^bits;

n = 8;

ef = SetPrecision[
   n p (betag  betah) (1 + ((n - 1)/16) (Q^2) (Sin[Pi/Q]^2)), prec];

ef2 = SetPrecision[
   n (p^2) (1/
      256) (((betag^2) (betah^2))) (512 (n + 
         1) + (32 (n - 1) (Q^2))/(Pi^2) + ((n - 
           1) (Q^2) (-32 Cos[(4 Pi)/Q] + 
           Pi (Sin[
               Pi/Q]^2) (16 (Pi + 4  n  Pi) + (n - 
                 2) Q ((n - 3) Pi Q (Sin[Pi/Q]^2) + 
                 16 Sin[(2 Pi)/Q]))))/(Pi^2)), prec];

a1 = SetPrecision[(ef2 - (ef^2)), prec];
a2 = SetPrecision[(ef2 - 5 (ef^2)), prec];
a3 = SetPrecision[-6 (ef^2), prec];

[Kappa] = SetPrecision[(-a2 + Sqrt[(a2^2) - 4 a1 a3])/(2 a1), prec];

[Theta] = SetPrecision[Sqrt[ef/([Kappa] ([Kappa] + 1))], prec];

f1 = SetPrecision[
   2^(-1 - [Kappa]/2) a b^(-[Kappa]/2) [Theta]^-[Kappa], prec];
f2 = SetPrecision[([Kappa] HypergeometricPFQ[{1/2 + [Kappa]/2, 
      1 + [Kappa]/2}, {3/2, 3/2 + [Kappa]/2}, 1/(2 b [Theta]^2)])/(
   Sqrt[2] Sqrt[b] [Theta] Gamma[(3 + [Kappa])/2]), prec];
f3 = SetPrecision[
   HypergeometricPFQ[{1/2 + [Kappa]/2, [Kappa]/2}, {1/2, 
     1 + [Kappa]/2}, 1/(2 b [Theta]^2)]/Gamma[1 + [Kappa]/2], prec];

SetPrecision[f1 (f3 - f2), prec]