I was looking to approximate this integral using ideas from Laplace's method, however I ended up with a $\color{red}{\text{divergent improper integral}}$. I want to know why such a method doesn't work in this case. I suspect it is because of how I define $k$ which might be interfering with the approximations.
For $\mu,\sigma>0$ and $m\in\mathbb{R}$ $$I=\int_0^\infty x^m\exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right)\ dx$$ my attempt is below.
The integrand is, $$ x^m\exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right)=\exp\left(m\log x-\frac{(x-\mu)^2}{2\sigma^2}\right) $$ which obtains its maximum value when the derivative of the exponent is $0$, $$\frac{m}{x}-\frac{x-\mu}{\sigma^2}=0\Rightarrow x=\frac{\mu+\sqrt{\mu^2+4\sigma m}}{2}$$ call the point $k$, then we have, $$\begin{align*} I&\sim\int_{k-\epsilon}^{k+\epsilon}\exp\left(m\log x-\frac{(x-\mu)^2}{2\sigma^2}\right)\ dx \\ &=\int_{-\epsilon}^\epsilon\exp\left(m\log(t+k)-\frac{(t+k-\mu)^2}{2\sigma^2}\right)\ dt \end{align*}$$ where we substituted $t=x-k$. By Taylor's formula, $$m\log(t+k)-\frac{(t+k-\mu)^2}{2\sigma^2}\sim m\log k-\frac{(k-\mu)^2}{2\sigma^2}+\left(\frac{m}{k}-\frac{\mu-k}{\sigma^2}\right)t+O(t^2)$$ when $t$ is in the neighborhood of $0$ (note that the first term is a constant). Keeping the first non-constant term, and exponentiating the rest via the Maclaurin series, $$\begin{align*}&\int_{-\epsilon}^\epsilon\exp\left(m\log k-\frac{(k-\mu)^2}{2\sigma^2}+\left(\frac{m}{k}-\frac{\mu-k}{\sigma^2}\right)t\right)(1+O(t^2))\ dt\\ \sim&\int_{-\epsilon}^\epsilon\exp\left(m\log k-\frac{(k-\mu)^2}{2\sigma^2}+\left(\frac{m}{k}-\frac{\mu-k}{\sigma^2}\right)t\right)\ dt \\ =&k^m\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\int_{-\epsilon}^\epsilon\exp\left(\left(\frac{m}{k}-\frac{\mu-k}{\sigma^2}\right)t\right)\ dt. \\ \end{align*}$$ Now changing the end points to $-\infty$ and $\infty$ only adds a subdominant function, $$k^m\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\color{red}{\int_{-\infty}^\infty\exp\left(\left(\frac{m}{k}-\frac{\mu-k}{\sigma^2}\right)t\right)\ dt}$$ however here, $$\frac{m}{k}-\frac{\mu-k}{\sigma^2}=\frac{m\sigma^2-uk+k^2}{k\sigma^2}=\frac{m\sigma^2+k^2}{k\sigma^2}-\frac{\mu}{\sigma^2}$$ with $$k:=\frac{\mu+\sqrt{\mu^2+4\sigma m}}{2}$$ seems to be positive for any $\mu,\sigma>0$ and $m\in\mathbb{R}$; this means the $\color{red}{\text{integral}}$ is divergent.
As @Gonçalo pointed out, the coefficient of $t^1$ was zero by the definition of $k$, so taking one more term in the expansion should resolve the issue.
By Taylor's formula,
$$\begin{align*}m\log(t+k)-\frac{(t+k-\mu)^2}{2\sigma^2}&\sim m\log k-\frac{(k-\mu)^2}{2\sigma^2}+\left(\frac{m}{k}-\frac{k-\mu}{\sigma^2}\right)t-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2+O(t^3) \\ &=m\log k-\frac{(k-\mu)^2}{2\sigma^2}-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2+O(t^3) \end{align*}$$ in the neighborhood of $t=0$. Keeping the first non-constant term and exponentiating the rest via the Maclaurin series, $$\begin{align*}I&\sim\int_{-\varepsilon}^\varepsilon \exp\left(m\log k-\frac{(k-\mu)^2}{2\sigma^2}-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2\right)(1+O(t^3))\ dt \\ &\sim k^m\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\int_{-\varepsilon}^\varepsilon\exp\left(-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2\right)\ dt \\ &\sim k^m\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\color{navy}{\int_{-\infty}^\infty\exp\left(-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2\right)\ dt} \\ &=2k^m\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\int_{0}^\infty\exp\left(-\left(\frac{m}{2k^2}+\frac{1}{2\sigma^2}\right)t^2\right)\ dt \end{align*}$$ note that the $\color{navy}{\text{integral}}$ is happily finite. Make a change of variables from $t$ to $s$ via, $$t=\dfrac{s}{\sqrt{\dfrac{m}{2k^2}+\dfrac{1}{2\sigma^2}}},\quad dt=\frac{ds}{\sqrt{\dfrac{m}{2k^2}+\dfrac{1}{2\sigma^2}}}$$ we find $$I\sim\frac{2k^m}{\sqrt{\dfrac{m}{2k^2}+\dfrac{1}{2\sigma^2}}}\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\int_{0}^\infty\exp(-s^2)\ ds$$ hence by the Gaussian integral, $$I\sim \frac{k^m\sqrt\pi}{\sqrt{\dfrac{m}{2k^2}+\dfrac{1}{2\sigma^2}}}\exp\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)$$ for large $k$.