But it is true:
You need a convenient locally convex space $E$ of functions in $x$ such that $g\mapsto \int_0^\infty g(x)\,dx$ is a bounded linear functional, and such that $y\mapsto f(\quad,y)$ is a differentiable curve in $E$, at least of Hoelder class $C^{1,\alpha}$. Then the interchange follows by the chain rule. The space $\mathcal S(-1,\infty)$ of smooth functions which decrease rapidly on the right suggests itself. This is a Frechet space with seminorms $sup \{(1+|x|^2)^m\partial_x^{k} g(x): -1<x<\infty\}$ for fixed $m,k$
Now we look at the paper: Peter W. Michor and David Mumford: A zoo of diffeomorphism groups on $\mathbb ℝ^n$. arXiv:1211.5704. (pdf)
By repeating the proof of 3.5 there (which uses the very strong theorem 2.4 of Frölicher and Kriegl) we see that:
The space $C^\infty(\mathbb R,\mathcal S(-1,\infty))$ of smooth curves in $\mathcal S(-1,\infty)$ consists of all functions $c\in C^\infty(\mathbb R\times(-1,\infty),\mathbb R)$ satisfying the following property:
$\bullet$ For all $k,m\in \mathbb N_{\ge0}$ and $\alpha\in \mathbb N_{\ge0}^n$, the expression $(1+|x|^2)^m\partial_t^{k}\partial^\alpha_xc(t,x)$ is uniformly bounded in $x\in (-1,\infty)$, locally in $t\in \mathbb R$.
Obviously, your function $f$ satisfies this.
EDIT: The "obviously" is wrong. $(1+x^2) y^3 e^{-y^2x}$ is not uniformly bounded for $(x,y)\in (-1,\infty)\times[-1,1]$. To see this insert the curve $x(t) = 1/(1+t^2)$, $y(t)=t$ to get $$(1+x^2) y^3 e^{-y^2x}= \Big(1+\frac{1}{(1+t^2)^2}\Big) t^3 e^{-\frac{t^2}{1+t^2}}, \quad t\in[-1.1]$$ which is unbounded.
Let me try again: $g\mapsto \int_0^\infty g(x)dx$ is a bounded linear functional on $L^1(-1,\infty)$, but $y\mapsto f(\quad,y)$ is not $Lip^1$ into (otherwise you could interchange).