Here comes 2016!

Alternate solution:

${ I }_{ n }=n\int _{ 0 }^{ 1 }{ { x }^{ n }{ e }^{ 2016x }dx }$

Now, we will upper-bound and lower-bound the integral. Since $x\le 1$ throughout the range of the integral, we have,

$I_n \leq n\int_{0}^{1} z^ne^{2016}dz\leq e^{2016}\hspace{20pt}(eqn1)$

Similarly, from integration by-parts we obtain $I_n=\frac{n}{n+1}e^{2016}-\frac{2016n}{n+1}\int_{0}^{1} e^{2016x}x^{n+1}dx\geq \frac{n}{n+1}e^{2016}-\frac{2016ne^{2016}}{(n+1)(n+2)}$

Taking limit as $n \to \infty$ , we have

$\lim_{n\to \infty}I_n\geq e^{2016}\hspace{40pt}(eqn2)$

The result now follows from $(eqn1)$ and $(eqn2)$ :

${ I }_{ n }=e^{2016}$

This solution bases on "Differentiation under integral sign" technique.

First $\displaystyle \lim _{ n\rightarrow \infty }{ \ln { \left( n\int _{ 0 }^{ 1 }{ \left( x+\frac { 2016{ x }^{ 2 } }{ n } \right) ^{ n }dx } \right) } } ={ \ln { \lim _{ n\rightarrow \infty }\left( n\int _{ 0 }^{ 1 }{ { x }^{ n }\left( 1+\frac { 2016{ x } }{ n } \right) ^{ n }dx } \right) } }$

Since $\displaystyle \lim _{ n\rightarrow \infty }\left( 1+\frac { 2016{ x } }{ n } \right) ^{ n }={ e }^{ 2016x }$

The integral becomes

$\displaystyle ={ \ln { \left( \lim _{ n\rightarrow \infty } n\int _{ 0 }^{ 1 }{ { x }^{ n }{ e }^{ 2016x }dx } \right) } }$

So, our problem now is to work out the above integral. Here is where Feynman's technique comes to play role.

Consider

$\displaystyle \int _{ 0 }^{ 1 }{ { e }^{ ax } } dx$

$\displaystyle =\frac { { e }^{ a }-1 }{ a }$

$\displaystyle \frac { { \partial }^{ n } }{ { \partial a }^{ n } } \int _{ 0 }^{ 1 }{ { e }^{ ax } } dx=\int _{ 0 }^{ 1 }{ { x }^{ n }{ e }^{ ax } } dx$

$\displaystyle \int _{ 0 }^{ 1 }{ { x }^{ n }{ e }^{ ax } } dx=\frac { { \partial }^{ n } }{ { \partial a }^{ n } } \left( \frac { { e }^{ a }-1 }{ a } \right)$

$\displaystyle \left( \frac { { e }^{ a }-1 }{ a } \right) =\sum _{ m=0 }^{ \infty }{ \frac { { a }^{ m } }{ \left( m+1 \right) ! } }$

Then

$\displaystyle \frac { { \partial }^{ n } }{ { \partial a }^{ n } } \left( \frac { { e }^{ a }-1 }{ a } \right) =\frac { { \partial }^{ n } }{ { \partial a }^{ n } } \sum _{ m=0 }^{ \infty }{ \frac { { a }^{ m } }{ \left( m+1 \right) ! } }$

$\displaystyle \frac { { \partial }^{ n } }{ { \partial a }^{ n } } \sum _{ m=0 }^{ \infty }{ \frac { { a }^{ m } }{ \left( m+1 \right) ! } } =\sum _{ m=0 }^{ \infty }{ \frac { { a }^{ m-n } }{ \left( m+1 \right) ! } } \frac { m! }{ \left( m-n \right) ! }$

Since every terms which has m less than n will be 0.Let's let $j=m-n$ and hence the sum will be reduced to

$\displaystyle \sum _{ m=0 }^{ \infty }{ \frac { { a }^{ m-n } }{ \left( m+1 \right) ! } } \frac { m! }{ \left( m-n \right) ! } =\sum _{ j=0 }^{ \infty }{ \frac { { a }^{ j } }{ \left( n+j+1 \right) j! } }$

$\displaystyle \lim _{ n\rightarrow \infty }{ n\sum _{ j=0 }^{ \infty }{ \frac { { a }^{ j } }{ \left( n+j+1 \right) j! } } } =\lim _{ n\rightarrow \infty }{ \sum _{ j=0 }^{ \infty }{ \frac { { na }^{ j } }{ \left( n+j+1 \right) j! } } } ={ \sum _{ j=0 }^{ \infty }{ \lim _{ n\rightarrow \infty } \frac { n }{ n+j+1 } \frac { { a }^{ j } }{ j! } } }$

$\displaystyle { \sum _{ j=0 }^{ \infty }{ \lim _{ n\rightarrow \infty } \frac { n }{ n+j+1 } \frac { { a }^{ j } }{ j! } } }={ \sum _{ j=0 }^{ \infty }{ \frac { { a }^{ j } }{ j! } } }={ e }^{ a }$

Substitute a=2016 as in the question we get

$\displaystyle \lim _{ n\rightarrow \infty }{ \ln { \left( n\int _{ 0 }^{ 1 }{ \left( x+\frac { 2016{ x }^{ 2 } }{ n } \right) ^{ n }dx } \right) } }=\ln{e^{2016}}=\boxed{2016}$

The function

$\displaystyle f(x) = \lim _{ n\rightarrow \infty }{ \left\{ n{ x }^{ n }\quad 0\le x\le 1\\ \quad \quad \quad \quad \quad \quad 0\quad elsewhere \right\} } = \delta{(x-1)}$

is a dirac delta function.

Properties of dirac delta function :

$\displaystyle \int _{ a-\epsilon }^{ a+\epsilon }{ \delta (x-a)dx } =1$ , for every non zero $\epsilon$ .

$\displaystyle \int _{ a-\epsilon }^{ a+\epsilon }{ \delta (x-a)f(x)dx } =f(a)$

So the integral here can be written as :

$\displaystyle \lim _{ n\rightarrow \infty }{ \int _{ 0 }^{ 1 }{ n{ x }^{ n }{ \left( 1+\frac { ax }{ n } \right) }^{ n }dx } } =\int _{ 0 }^{ 1 }{ \delta (x-1){ e }^{ ax }dx } ={ e }^{ a }$

Hence we are done put $\displaystyle a=2016$ to get the answer :