We have the joint PDF for $X$ and $Y$ as: $$\frac{1}{2\pi} \cdot \frac{1}{\sqrt{a^2b^2-c^2}} \exp\left\{\frac{-(a^2(x-m)^2 - 2c(x-m)(y-n)+b^2(y-n)^2)}{2(a^2b^2 - c^2)}\right\} dx dy$$
Where $m,n \in\mathbb{R}, a>0, b>0, c\in\mathbb{R}, |c|<ab$.
Now I am trying to find Cov($X,Y$), only through integration of transcendental functions or through properties of Variance, essentially at an undergrad level of Calculus.
The way I've been attempting to do it is by replacing variables, completing the square, and getting:
$$\frac{1}{2\pi\sqrt d a ^2b^2} \cdot \iint_{\mathbb{R^2}}\exp\left\{\frac{-(\frac{t^2}{2ab}+\frac{w^2}{2d})}{b}\right\}\cdot {\left\{tw+\frac{ct^2}{ab}\right\}} dt dw$$
by making the following substitutions:
$d=a^2b^2-c^2, (x-m) = u, (y-n) = v, au = z, bv = t$, and finally, $w = z - \frac{ct}{ab}$.
Now applying all of these transformations, I get the double integral above. But now I am confused as to what to do. Also, I may be overthinking this question, and making there is an easier way through polar coordinates that I can reach the same result? If you had any feedback it would be appreciated. Thank you for your time.
You seem to have some algebra mistakes in your calculation, leading to a wrong answer. A cleaner set of substitutions is: $$ z:=\frac{x-m}b,\quad t:=\frac{y-n}a,\quad\rho:=\frac c{ab}.\tag1 $$ Assuming you have established that $E(X)=m$ and $E(Y)=n$, the covariance between $X$ and $Y$ is $$ \operatorname{Cov}(X,Y)=\iint (x-m)(y-n)f(x,y)\,dxdy\tag2. $$ Applying the substitutions (1) you will get $$ \begin{align} &\iint bz\, at\, f(bz+m,at+n)\,bdz\, a dt\\ & = (ab)^2\iint zt \frac1{2\pi ab\sqrt{1-\rho^2}}\exp\left\{-\frac{a^2b^2z^2-2cabzt+b^2a^2t^2}{2(a^2b^2-c^2)}\right\}\,dzdt\\ & = ab\iint zt \frac1{2\pi \sqrt{1-\rho^2}}\exp\left\{-\frac{z^2-2\rho zt+t^2}{2(1-\rho^2)}\right\}\,dzdt\\ & = ab\iint zt\frac1{\sqrt{2\pi(1-\rho^2)}}\exp\left\{-\frac{(z-\rho t)^2}{2(1-\rho^2)}\right\}\frac1{\sqrt{2\pi}}\exp\left\{-\frac{t^2}2\right\}\,dzdt.\tag3 \end{align} $$ To evaluate (3), use your substitution $w:=z-\rho t$ to obtain $$ ab\iint (w+\rho t)t\frac1{\sqrt{2\pi(1-\rho^2)}}\exp\left\{-\frac{w^2}{2(1-\rho^2)}\right\}\frac1{\sqrt{2\pi}}\exp\left\{-\frac{t^2}2\right\}\,dwdt.\tag4 $$ Breaking (4) into two pieces, the first term is $$ ab\left(\int w\frac1{\sqrt{2\pi(1-\rho^2)}}\exp\left\{-\frac{w^2}{2(1-\rho^2)}\right\}\,dw\right)\left( \int t\frac1{\sqrt{2\pi}}\exp\left\{-\frac{t^2}2\right\}\,dt\right)\tag5 $$ which equals $ab\cdot 0\cdot0$, and the second term is $$ \rho ab\left(\int \frac1{\sqrt{2\pi(1-\rho^2)}}\exp\left\{-\frac{w^2}{2(1-\rho^2)}\right\}\,dw\right)\left( \int t^2\frac1{\sqrt{2\pi}}\exp\left\{-\frac{t^2}2\right\}\,dt\right)\tag6 $$ which equals $\rho ab\cdot1\cdot1=c$.