Laplace method (or other integral asymptotic) with near-corner Consider the integral
$$\int_{-\infty}^\infty \exp(-\sqrt{h^2+M^2x^2}) dx.$$
Here $h$ is a small positive parameter and $M$ is a large positive parameter. I would like to obtain a "reasonably uniform" asymptotic approximation for this integral in the limit of large $M$ and small $h$, specifically when $h$ goes to zero before $M$ goes to infinity. 
The difficulty is that the leading order part of the Laplace method sees $\sqrt{h^2+M^2 x^2}$ as $h+\frac{M^2}{2h} x^2$, a quadratic function, but in fact this approximation is only any good where $|x| \ll h/M$. By contrast there is a significant contribution to the integration over an interval of length on the order of $1/M$, which is much larger. Higher order Taylor approximations never see this because they just keep on assuming that $|x| \ll h/M$ and thus proceed to divide by larger and larger powers of $h$. 
An obvious alternative is to sacrifice accuracy in this $O(h/M)$ vicinity of $0$, for example by suppressing $h^2$ altogether, but this obviously does not achieve $o(h)$ accuracy, which is required for my application. Is there another workaround for this situation? Perhaps by "matching" the two approximations which are valid in different regimes?

As noticed, a simple Laplace method cannot be used here, as 2 scales are involved. A uniform asymptotic expansion should be found. Alternatively, in this case, we can recognize a modified Bessel function. Indeed, changing $x=\frac{h}{M}\sinh t$, the integral can be written as
\begin{align}
 I&=\int_{-\infty}^\infty \exp(-\sqrt{h^2+M^2x^2}) \,dx\\
 &=\frac{h}{M}\int_{-\infty}^\infty \exp(-h\cosh t)\cosh t \,dt
\end{align}
which is proportional to an integral representation of a modified Bessel function (DLMF):
\begin{equation}
 K_{\nu}\left(z\right)=\int_{0}^{\infty}e^{-z\cosh t}\cosh\left(\nu t\right)%
\mathrm{d}t
\end{equation} 
with $\nu=1,z=h$,
\begin{equation}
 I=\frac{2h}{M}K_1(h)
\end{equation} 
Using the series expansion of the Bessel function near $h=0$ (DLMF),
\begin{align}
 I&\sim \frac{2h}{M}\left[ \frac{1}{h}+\frac{h}{2}\left(\ln\frac{h}{2} +\gamma-\frac{1}{2}\right)+\ldots\right]\\
 &\sim\frac{2}{M}+\frac{h^2}{M}\left(\ln\frac{h}{2} +\gamma-\frac{1}{2}\right)+\ldots
\end{align} 
EDIT:Another method using the Mellin transform technique
Changing $x=th/M$ and then $u=\sqrt{1+t^2}-1$, the problem is equivalent to find the small $h$ behavior of 
\begin{align}
 I&=2\frac{h}{M}\int_0^\infty\exp(-h\sqrt{1+t^2})\,dt\\
 &=2\frac{h}{M}e^{-h}\int_0^\infty\frac{u+1}{\sqrt{u(u+2)}}\exp(-hu)\,du
\end{align}
We have thus to find a Laplace transforms with a small parameter. A classical method which uses the Mellin transform technique is given in (DLMF). Intermediate results are given below, with the help of a CAS. Defining
\begin{equation}
 H(u)=\frac{u+1}{\sqrt{u(u+2)}}
\end{equation} he following behaviors hold:\
 \begin{array}{lll}
 H(u)&\sim 1+\frac{1}{2u^2}+O\left( u^{-3} \right)& \text{ for }u\to\infty\\
 &\sim O(u^{-1/2})&  \text{ for } u\to 0
\end{array}
and the Mellin transform is
\begin{equation}
 \mathcal{M}\left[H(u) \right](z)=\pi^{-1/2}2^{z-1}(z-1)\Gamma(-z)\Gamma\left( z-\frac{1}{2} \right)
\end{equation} 
For the function 
\begin{equation}
 F(z)=-h^{z-1}\Gamma(1-z)\mathcal{M}\left[H(u) \right](z)
\end{equation} 
the residues for $z=0,1,2,3$ are
 \begin{align}
 \left. \operatorname{res}F(z)\right|_{z=0}&=\frac{1}{h}\\
\left. \operatorname{res}F(z)\right|_{z=1}&=1\\
 \left. \operatorname{res}F(z)\right|_{z=2}&=\frac{h}{2}\left[ \ln\frac{h}{2}+\gamma+\frac{1}{2}\right]\\
 \left. \operatorname{res}F(z)\right|_{z=3}&=\frac{h^2}{12}\left[6 \ln\frac{h}{2}+6\gamma-1\right]
 \end{align}
With $e^{-h}= 1-h+h^2/2+O(h^3)$,
\begin{align}
 I&\sim \frac{2h}{M}(1-h+\frac{h^2}{2})\left[ \frac{1}{h}+\frac{h}{2}\left( \ln\frac{h}{2}+\gamma+\frac{1}{2}\right)+\frac{h^2}{12}\left(6 \ln\frac{h}{2}+6\gamma-1\right)\right]\\
 &\sim \frac{2}{M}+\frac{h^2}{M}\left( \ln\frac{h}{2}+\gamma-\frac{1}{2} \right)
\end{align}
