I am looking for a reference (or a somewhat simple proof) for the following result, which for instance Mathematica spits out without too much effort. Here $a,b,c \in \mathbb{R}$ are constants satisfying $a, c < 0$ and $b^2 < 4 a c$.
$$\int_0^{\infty}\int_0^{\infty} \exp(a x^2 + b x y + c y^2) \, dx \, dy = \frac{1}{2\sqrt{4ac-b^2}} \left(\pi + 2 \arctan\left(\frac{b}{\sqrt{4ac-b^2}}\right)\right).$$
One idea I had was to write:
\begin{align} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \exp(a x^2 + b x y + c y^2) \, dx \, dy &= 2 \int_0^{\infty}\int_0^{\infty} \exp(a x^2 + b x y + c y^2) \, dx \, dy \\ &+ 2 \int_0^{\infty}\int_0^{\infty} \exp(a x^2 - b x y + c y^2) \, dx \, dy \end{align}
The integral on the left can easily be computed as it corresponds to the total probability of the corresponding bivariate Gaussian distribution (minus scaling factors), but then I'm stuck at solving the other integral with $(-b)$ substituted for $b$, which seems just as complicated to compute and is not easily relatable to the original integral. So I'm not sure if this approach leads anywhere.
A more direct approach would be to simply compute both integrals explicitly, one at a time. Doing one integral first, a human can verify without too much effort that
$$\int_0^{\infty} \exp(a x^2 + b x y + c y^2) \, dx = \sqrt{\frac{\pi}{-4 a}} \cdot \exp\left(\frac{(4 a c -b^2) y^2}{a}\right) \left(\text{erf}\left(\frac{b y}{2 \sqrt{-a}}\right)+1\right).$$
Doing the second integral with $\exp(u y^2) (\text{erf}(v y) + 1)$ is not so straightforward though, and for instance I could not find a reference for computing such integrals (and getting an arctangent in the process) in Abramowitz and Stegun's handbook. So if someone has a reference for integrating $\exp(u y^2) \text{erf}(v y)$ over $y > 0$ for $u, v$ as above, that would also be appreciated.