fxy = Exp[-1*(((x - ux)/sx)^2 + ((y - uy)/sy)^2 -
2 r ((x - ux)*(y - uy)/(sx*sy)))/(2 (1 - r^2))]/(2*Pi*sx*
sy*(1 - r^2)^0.5);
fy = (1/(sy (2*Pi)^0.5))*Exp[-(1/2)*((y - uy)/sy)^2];
f = Rationalize[x^2*y*(fxy/fy), 0] // Simplify // PowerExpand
(*(35567099 E^((sy (ux - x) + r sx (-uy + y))^2/(
2 (-1 + r^2) sx^2 sy^2)) x^2 y)/(89153496 Sqrt[1 - r^2] sx)*)
INT1 = Integrate[f, {y, -Infinity, Infinity},
Assumptions -> {ux > 0, uy > 0, sx > 0, sy > 0, 0 < r < 1,
x \[Element] Reals}]
(*(35567099 Sqrt[\[Pi]/2] sy x^2 (r sx uy +
sy (-ux + x)))/(44576748 r^2 sx^2)*)
Integrate[INT1, {x, -Infinity, Infinity},
Assumptions -> {ux > 0, uy > 0, sx > 0, sy > 0, 0 < r < 1}]
(*"Integral of 1 does not converge on {-\[Infinity],\[Infinity]}"*)
Looks like integral is divergent.