Use Laplace's method with $\int_{0}^{\infty} e^{x(3u-u^3)}du$ as $x\rightarrow \infty$ Use Laplace's method with $\int_{0}^{\infty} e^{x(3u-u^3)}du$ as $x\rightarrow \infty$.
I'm confused about how to taylor expand about u=1? How do I continue? Obviously first of all I have converted it to:
$$\int_{1-\epsilon}^{1+\epsilon}e^{x(3u-u^3)}du$$
but what now?

Let's see the where the usual reasoning for Laplace's method leads us.
Let me write the integral as $\exp{(x f(u))}$, where $f(u)$ has a global maximum at $u_* = 1$. We could anticipate that the major contribution to the integral is coming from the neighbourhood of $u = u_*$ as $x \to \infty$. If we expanded the argument of the exponential about $u = u_*$ we would find:
$$ I(x) = \int^\infty_0 \mathrm{e}^{xf(u)} \, \mathrm{d}u \approx \int^\infty_0 \mathrm{e}^{x \left[ f(u_*) + \frac{1}{2}f''(u_*)(u-u_*)^2 \right] } \, \mathrm{d}u, $$ where $f(u_*) = 2$ and $f''(u_*) = -6$ (of course, $f'(u_*)=0$). Now, we call Mathematica to solve the last integral to come up with:
$$I(x) \sim \sqrt{\frac{\pi }{12 x}} \mathrm{e}^{2 x} \left(1+\text{erf}\left( \sqrt{3x}\right)\right), \quad x \to \infty $$ where $\mathrm{erf}$ is the error function. We have obtained two terms in the asymptotic expansion of $I(x)$. 
Hope this helps

Note that if you do not retain the second order in the expansion of $f$ you end up with a non-convergent integral for $u$ leading to the value of $I$ at $x \to \infty$, which is not precisely what we desire.

Here's a comparison between the numerical integration of $I(x)$ (black solid line) and the two-term asymptotic expansion (dashed line):

