<s>[INST] Double Integration Problem $\int_{0}^{1} \int_0^1 \frac{1}{1+y(x^2-x)}dydx$ 
Compute
  $$I = \int_{0}^{1} \int_0^1 \frac{1}{1+y(x^2-x)}dydx$$

Here are my steps:
$$\begin{split}
I &=\int_{0}^{1} \left(\int_0^1 \frac{dy}{1+y(x^2-x)}\right)dx\\
  &=\int_{0}^{1} \left[\frac{\ln(1+y(x^2-x))}{x^2-x}\right]_0^1dx\\
  &=\int_{0}^{1} \left[\frac{\ln(1+(1)(x^2-x))}{x^2-x}
                      -\frac{\ln(1+(0)(x^2-x))}{x^2-x}\right]dx\\
  &=\int_{0}^{1} \frac{\ln(1+x^2-x)}{x^2-x}dx
\end{split}
$$
And here I can't find any substitution to solve this integral.
Can anyone help me?
By the way, I also used Simpson's 3/8 method to find the approximation and got $1.063$.
But I want to find it using Calculus.
 [/INST] Here is an alternate way that avoids your integral.
Writing the integrand in terms of a geometric sum we have:
$$\frac{1}{1 + y(x^2 - x)} = \sum_{n = 0}^\infty y^n (x - x^2)^n, \quad |x|, |y| < 1,$$
then
\begin{align}
I &= \int_0^1 \int_0^1 \frac{1}{1 + y(x^2 - x)} \, dy dx\\
&= \sum_{n = 0}^\infty \int_0^1 y^n \, dy \int_0^1 x^n (1 - x)^n \, dx\\
&= \sum_{n = 0}^\infty \frac{1}{n + 1} \cdot \operatorname{B}(n + 1, n + 1)\\
&= \sum_{n = 0}^\infty \frac{1}{n + 1} \cdot \frac{(n!)^2}{(2n + 1)!}\\
&= \sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}}.\tag1
\end{align}
Here $\operatorname{B}(x,x)$ denoted the Beta function while $\binom{2n}{n} = \frac{(2n)!}{(n!)^2}$ is the central binomial coefficient.
To find the sum in (1), the following (reasonably?) well known Maclaurin series can be recalled (for a proof, see here)
$$\frac{\sin^{-1} x}{\sqrt{1 - x^2}} = \sum_{n = 0}^\infty \frac{2^{2n} x^{2n + 1}}{(2n + 1) \binom{2n}{n}}, \quad |x| < 1.$$
Enforcing a substitution of $x \mapsto \sqrt{x}$, after rearranging one has
$$\frac{\sin^{-1} \sqrt{x}}{\sqrt{x} \sqrt{1 - x}} = \sum_{n = 0}^\infty \frac{2^{2n} x^n}{(2n + 1) \binom{2n}{n}}, \quad |x| < 1.$$
Integrating both sides with respect to $x$ from $0$ to $\frac{1}{4}$ leads to
$$\sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}} = 4 \int_0^{\frac{1}{4}} \frac{\sin^{-1} \sqrt{x}}{\sqrt{x} \sqrt{1 - x}} \, dx.$$
Making a substitution of $x = \sin^2 u$ in the integral readily leads to
\begin{align}
\sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}} &= 8 \int^{\frac{\pi}{6}}_0 u \, du = \frac{\pi^2}{9},
\end{align}
from which it follows that:
$$\int_0^1 \int_0^1 \frac{1}{1 + y(x^2 - x)} \, dy dx = \frac{\pi^2}{9}.$$
</s>