Theorem: First Order ODE/dx = (y over (1 - x^2 y^2)) dx + (x over (1 - x^2 y^2)) dy

is an exact differential equation with solution:
:$\map \ln {\dfrac {1 + x y} {1 - x y} } - 2 x = C$

Proof:
First express $(1)$ in the form:
:$(2): \quad \paren {\dfrac y {1 - x^2 y^2} - 1} + \paren {\dfrac x {1 - x^2 y^2} } \dfrac {\d y} {\d x}$
Let:
:$\map M {x, y} = \dfrac y {1 - x^2 y^2} - 1$
:$\map N {x, y} = \dfrac x {1 - x^2 y^2}$
Then:
{{begin-eqn}}
{{eqn | l = \frac {\partial M} {\partial y}
      | r = \frac {\paren {1 - x^2 y^2} - y \paren {- 2 x^2 y} } {\paren {1 - x^2 y^2}^2}
      | c = Quotient Rule for Derivatives
}}
{{eqn | r = \frac {1 + x^2 y^2} {\paren {1 - x^2 y^2}^2}
      | c = 
}}
{{eqn | l = \frac {\partial N} {\partial x}
      | r = \frac {1 + x^2 y^2} {\paren {1 - x^2 y^2}^2}
      | c = similarly
}}
{{end-eqn}}
Thus $\dfrac {\partial M} {\partial y} = \dfrac {\partial N} {\partial x}$ and $(1)$ is seen by definition to be exact.
By Solution to Exact Differential Equation, the solution to $(1)$ is:
:$\map f {x, y} = C$
where:
{{begin-eqn}}
{{eqn | l = \dfrac {\partial f} {\partial x}
      | r = \map M {x, y}
      | c = 
}}
{{eqn | l = \dfrac {\partial f} {\partial y}
      | r = \map N {x, y}
      | c = 
}}
{{end-eqn}}
Hence:
{{begin-eqn}}
{{eqn | l = f
      | r = \map M {x, y} \rd x + \map g y
      | c = 
}}
{{eqn | r = \int \paren {\dfrac y {1 - x^2 y^2} - 1} \rd x + \map g y
      | c = 
}}
{{eqn | r = \int \dfrac {y \rd x} {1 - x^2 y^2} - \int \rd x + \map g y
      | c = Integral Operator is Linear
}}
{{end-eqn}}
Substitute $z = x y$ to obtain:
:$\dfrac {\d z} {\d x} = y$
This gives:
{{begin-eqn}}
{{eqn | l = \int \frac {y \rd x} {1 - x^2 y^2} - \int \rd x + \map g y
      | r = \int \frac {\rd z} {1 - z^2} - \int \rd x + \map g y
      | c = substituting $z = x y$
}}
{{eqn | r = \int \frac 1 2 \paren {\frac 1 {1 + z} + \frac 1 {1 - z} } \rd z - \int \rd x + \map g y
      | c = Integration by Partial Fractions
}}
{{eqn | r = \frac 1 2 \int \dfrac {\rd z} {1 + z} + \frac 1 2 \int \dfrac {\rd z} {1 - z} - \int \rd x + \map g y
      | c = Integral Operator is Linear
}}
{{eqn | r = \frac 1 2 \map \ln { 1 + z } - \frac 1 2 \map \ln { 1 - z } - \int \rd x + \map g y
      | c =
}}
{{eqn | r = \frac 1 2 \map \ln {\frac {1 + z} {1 - z} } - \int \rd x + \map g y
      | c = Difference of Logarithms
}}
{{eqn | ll= \leadsto
      | l = \map M {x, y} \rd x + \map g y
      | r = \frac 1 2 \map \ln {\frac {1 + x y} {1 - x y} } - x + \map g y
      | c = 
}}
{{end-eqn}}
and:
{{begin-eqn}}
{{eqn | l = f
      | r = \int \map N {x, y} \rd y + \map h x
      | c = 
}}
{{eqn | r = \int \frac {x \rd y} {\paren {1 - x^2 y^2} }  + \map h x
      | c = 
}}
{{eqn | r = \int \frac {\rd z} {\paren {1 - z^2} } + \map h x
      | c = substituting $z = x y$
}}
{{eqn | r = \int \frac 1 2 \paren {\frac 1 {1 + z} + \frac 1 {1 - z} } \rd z + \map h x
      | c = Integration by Partial Fractions
}}
{{eqn | r = \frac 1 2 \int \dfrac {\rd z} {1 + z} + \frac 1 2 \int \dfrac {\rd z} {1 - z} + \map h x
      | c = Integral Operator is Linear
}}
{{eqn | r = \frac 1 2 \map \ln { 1 + z } - \frac 1 2 \map \ln { 1 - z } + \map h x
      | c =
}}
{{eqn | r = \frac 1 2 \map \ln {\frac {1 + z} {1 - z} } + \map h x
      | c = Difference of Logarithms
}}
{{eqn | r = \frac 1 2 \map \ln {\frac {1 + x y} {1 - x y} } + \map h x
      | c =
}}
{{end-eqn}}
Thus:
:$\map f {x, y} = \dfrac 1 2 \map \ln {\dfrac {1 + x y} {1 - x y} } - x$
and by Solution to Exact Differential Equation, the solution to $(1)$, after simplification, is:
:$\map \ln {\dfrac {1 + x y} {1 - x y} } - 2 x = C$
Note that from the {{Defof|Natural Logarithm}}, $\map \ln z$ is only defined when $z > 0$.
This means $\dfrac {1 + x y} {1 - x y} > 0$
First assume $x y > 0$:
:$\ds \dfrac {1 + x y} {1 - x y} > 0 \to \paren {1 + x y} > 0 \to x y > -1$
Next assume $x y < 0$:
:$\ds \dfrac {1 + -\paren {x y} } {1 - -\paren {x y} } > 0 \to \paren {1 - x y} > 0 \to x y < 1$
Therefore, $\size {x y} < 1$
{{qed}}