Answer
We evaluate $I = \mathop \smallint \limits_1^3 \mathop \smallint \limits_0^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x$ using two methods.
Method 1. Hold $x$ constant and evaluate the inner integral with respect to $y$
Method 2. Change the order of integration: hold $y$ constant and integrate first with respect to $x$
The results are the same:
$I = \mathop \smallint \limits_1^3 \mathop \smallint \limits_0^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x \simeq 4.64356$.
Using Fubini's Theorem to change the order of integration (Method 2) is easier.
Work Step by Step
Evaluate $I = \mathop \smallint \limits_1^3 \mathop \smallint \limits_0^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x$.
Method 1. Hold $x$ constant and evaluate the inner integral with respect to $y$
$I = \mathop \smallint \limits_{x = 1}^3 \left( {\mathop \smallint \limits_{y = 0}^1 y{{\rm{e}}^{xy}}{\rm{d}}y} \right){\rm{d}}x$
Write the inner integral:
$S\left( x \right) = \mathop \smallint \limits_{y = 0}^1 y{{\rm{e}}^{xy}}{\rm{d}}y$
Let $u=y$ and $v = {{\rm{e}}^{xy}}$. So, $du = dy$ and $dv = x{{\rm{e}}^{xy}}dy$.
Thus,
$S\left( x \right) = \frac{1}{x}\mathop \smallint \limits_{y = 0}^1 u{\rm{d}}v$
Recall the Integration by Parts of Section 8.1,
$\smallint u{\rm{d}}v = uv - \smallint v{\rm{d}}u$
Using the Integration by Parts, we obtain
$S\left( x \right) = \frac{1}{x}\left( {uv|_{y = 0}^1 - \mathop \smallint \limits_{y = 0}^1 v{\rm{d}}u} \right)$
$ = \frac{1}{x}\left( {y{{\rm{e}}^{xy}}|_{y = 0}^1 - \mathop \smallint \limits_{y = 0}^1 {{\rm{e}}^{xy}}{\rm{d}}y} \right)$
$ = \frac{1}{x}\left( {{{\rm{e}}^x} - \frac{1}{x}{{\rm{e}}^{xy}}|_{y = 0}^1} \right)$
$ = \frac{1}{x}\left( {{{\rm{e}}^x} - \frac{1}{x}\left( {{{\rm{e}}^x} - 1} \right)} \right)$
$ = {{\rm{e}}^x}{x^{ - 1}} - {{\rm{e}}^x}{x^{ - 2}} + {x^{ - 2}}$
So,
$I = \mathop \smallint \limits_{x = 1}^3 \left( {\mathop \smallint \limits_{y = 0}^1 y{{\rm{e}}^{xy}}{\rm{d}}y} \right){\rm{d}}x$
$ = \mathop \smallint \limits_{x = 1}^3 \left( {{{\rm{e}}^x}{x^{ - 1}} - {{\rm{e}}^x}{x^{ - 2}} + {x^{ - 2}}} \right){\rm{d}}x$
$ = \mathop \smallint \limits_{x = 1}^3 \left( {{{\rm{e}}^x}{x^{ - 1}} - {{\rm{e}}^x}{x^{ - 2}}} \right){\rm{d}}x + \mathop \smallint \limits_1^3 {x^{ - 2}}{\rm{d}}x$
Using the given formula: $\smallint {{\rm{e}}^x}\left( {{x^{ - 1}} - {x^{ - 2}}} \right){\rm{d}}x = {x^{ - 1}}{{\rm{e}}^x} + C$ yields
$I = {x^{ - 1}}{{\rm{e}}^x}|_1^3 - {x^{ - 1}}|_1^3$
$ = \frac{1}{3}{{\rm{e}}^3} - {\rm{e}} - \frac{1}{3} + 1$
$ = \frac{1}{3}{{\rm{e}}^3} - {\rm{e}} + \frac{2}{3}$
$ \simeq 4.64356$
So, $I = \mathop \smallint \limits_1^3 \mathop \smallint \limits_0^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x \simeq 4.64356$.
Method 2. Change the order of integration: hold $y$ constant and integrate first with respect to $x$
Using Fubini's Theorem we change the order of integration and re-write the integral $I = \mathop \smallint \limits_{x = 1}^3 \mathop \smallint \limits_{y = 0}^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x$
$I = \mathop \smallint \limits_{y = 0}^1 \left( {\mathop \smallint \limits_{x = 1}^3 y{{\rm{e}}^{xy}}{\rm{d}}x} \right){\rm{d}}y$
$ = \mathop \smallint \limits_{y = 0}^1 \left( {{{\rm{e}}^{xy}}|_1^3} \right){\rm{d}}y$
$ = \mathop \smallint \limits_{y = 0}^1 \left( {{{\rm{e}}^{3y}} - {{\rm{e}}^y}} \right){\rm{d}}y$
$ = \left( {\frac{1}{3}{{\rm{e}}^{3y}} - {{\rm{e}}^y}} \right)|_0^1$
$ = \frac{1}{3}{{\rm{e}}^3} - {\rm{e}} - \frac{1}{3} + 1$
$ = \frac{1}{3}{{\rm{e}}^3} - {\rm{e}} + \frac{2}{3}$
$ \simeq 4.64356$
So, $I = \mathop \smallint \limits_1^3 \mathop \smallint \limits_0^1 y{{\rm{e}}^{xy}}{\rm{d}}y{\rm{d}}x \simeq 4.64356$.
Using Fubini's Theorem to change the order of integration method is easier.