The error function $\mathrm{erf}(x)$ is defined as: $$\mathrm{erf}(x):=\frac {2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^{2}}dt\tag{1}$$
Let us define the following 3 functions:
$$h(y)=\frac{1}{2}\left(\mathrm{erf}(1+i y)+\mathrm{erf}(1-i y)\right)\tag{2}$$
$$f(y)=e\sqrt{\pi}(2y^2+1)e^{-y^2}h(y)+2\left(\cos(2y) - y\sin(2y)\right) \tag{3}$$
$$g(y)=e\sqrt{\pi}y(2y^2-1)e^{-y^2}h(y)+2\left(y\cos(2y)+(y^2-1)\sin(2y)\right) \tag{4}$$
We also have: $$\frac{df}{dy}(y)=-2 g(y)$$
Here is a plot of $\color{red}{\log(f^2(y))}$ and $\color{blue}{\log(g^2(y))}$ vs. $y$:

Question: How to prove that the following sum-of-square function $\phi(y)$ is always positive:
$$\phi(y):=f^2(y)+g^2(y)>0,\qquad y\in \mathbb{R}\tag{5}$$
This problem showed up when I was looking for some approximation to Riemann $\Xi(z)$ function.
I will soon describe several unsuccessful attempts that I tried.
Best regards- mike