Solution:

We have,

$\begin{aligned} I &=\int_{0}^{1} \int_{x^{2}}^{x} x y(x+y) d y \ d x \\ I &=\int_{0}^{1}\left[\frac{x^{2} y^{2}}{2}+\frac{x y^{3}}{3}\right]_{x^{2}}^{x} d x \end{aligned}$

$\begin{aligned} I &=\int_{0}^{1}\left[\frac{x^{4}}{2}+\frac{x^{4}}{3}-\frac{x^{6}}{2}-\frac{x^{7}}{3}\right] d x \\ I &=\left[\frac{5}{6} \cdot \frac{x^{5}}{5}-\frac{x^{7}}{14}-\frac{x^{8}}{24}\right]_{0}^{1} \end{aligned}$

$\begin{aligned} I &=\frac{1}{6}-\frac{1}{14}-\frac{1}{24} \\ I &=\frac{3}{56} \end{aligned}$