1.) Yes, but the result is false. Correctly:
$$
\sqrt{\frac{2\pi}M} \cdot e^{M g(x_0)} \approx \int \sqrt{-g''(x)}\cdot e^{Mg(x)}\,\mathrm{d}x
$$
Which is something, but it is not trivial for me how an approximation for ${g(x_0)}$ follows. You have to be careful with the estimations.
2.) 3.) 4.) I think that it doesn't work for those cases. Here is why.
The reason the Laplace's method works is that $\lvert e^x\rvert$ is very small when $\mathbb{R}\ni x\to-\infty$. If $g(x)$ has a global maximum at $x_0$ then $e^{M g(x_0)}\cdot \int \frac{e^{M g(x)}}{e^{M g(x_0)}}$ uses the property that the integrand is small everywhere except near the global maximum.
So the real, global maximum is required (however multivariate generalizations exist).
Note that in practice, it is harder to evaluate the integral. On that I advise statistical physics papers, for example: http://www3.eng.cam.ac.uk/research_db/publications/gc121
At the end you might find that if you have enough points to approximate the integral then you have enough points to take just the maximum. So the Laplace's method is superfluous.
$$I=\underbrace{\int_0^{T}\frac{d\tau}{\sqrt{T^2-\tau^2}}\exp\left(-\frac{a^2}{4(T+\tau)}-\frac{b^2}{4(T-\tau)}\right)}\\{\tau\rightarrow \frac{T(1-y^2)}{1+y^2}}$$
$$=\exp\left(-\frac{a^2+b^2}{8T}\right)\int_0^{1} \frac{dy}{1+y^2}\exp\left(-\frac{a^2y^2+\frac{b^2}{y^2}}{8T}\right)$$
Let's parametrize the integral in order to build two ODE: $$ f(\alpha,\beta)=\int_0^{1} \frac{dy}{1+y^2}\exp\left(-\alpha^2y^2-\frac{\beta^2}{y^2}\right)$$
$$g(\alpha,\beta)=f(\alpha ,\beta)-\frac{1}{2\alpha}\frac{\partial f}{\partial\alpha}=\int_{0}^{1}\exp\left(-\alpha^2y^2-\frac{\beta^2}{y^2}\right)dy$$
Now, let's use $g(\alpha,\beta)$ to solve the remaining integral: $$\alpha g(\alpha,\beta)+\frac{1}{2}\frac{\partial g}{\partial\beta}=\exp\left(2\alpha\beta\right)\int_{0}^{1}\left(\alpha-\frac{\beta}{y^2}\right)\exp\left[-\left(\alpha y+\frac{\beta}{y}\right)^2 \right]dy=-\frac{\sqrt{\pi}}{2}\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)$$
$$\frac{\partial g}{\partial\beta}+2\alpha g(\alpha,\beta)=-\sqrt{\pi}\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)$$
$$\mu_1(\beta)=\exp\left(\int 2\alpha d\beta\right)=\exp\left(2\alpha\beta\right)$$
$$\exp\left(2\alpha\beta\right)\frac{\partial g}{\partial\beta}+2\alpha\exp\left(2\alpha\beta\right) g(\alpha,\beta)=-\sqrt{\pi}\exp\left(4\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)$$
$$\left(\exp\left(2\alpha\beta\right)g(\alpha,\beta)\right)'=-\sqrt{\pi}\int \exp\left(4\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)d\beta=-\frac{\sqrt{\pi}}{4\alpha}\left(\exp\left(4\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\operatorname{erf}\left(\alpha-\beta\right)\right)$$
$$g\left(\alpha,\beta\right)=-\frac{\sqrt{\pi}}{4\alpha}\left(\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\exp\left(-2\alpha\beta\right)\ \operatorname{erf}\left(\alpha-\beta\right)\right)+c$$
$$g\left(\alpha,0\right)=\int_{0}^{1}\exp\left[-\left(\alpha y\right)^2 \right]dy=\frac{\sqrt{\pi}}{2\alpha}\ \operatorname{erf}\left(\alpha\right)=-\frac{\sqrt{\pi}}{4\alpha}\left(\ \operatorname{erfc}\left(\alpha\right)-\ \operatorname{erf}\left(\alpha\right)\right)+c\\ c=\frac{\sqrt{\pi}}{4\alpha}$$
$$\boxed{g\left(\alpha,\beta\right)=-\frac{\sqrt{\pi}}{4\alpha}\left(\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\exp\left(-2\alpha\beta\right)\ \operatorname{erf}\left(\alpha-\beta\right)-1\right)}$$
Using this result in the first ODE: $$f(\alpha ,\beta)-\frac{1}{2\alpha}\frac{\partial f}{\partial\alpha}=-\frac{\sqrt{\pi}}{4\alpha}\left(\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\exp\left(-2\alpha\beta\right)\ \operatorname{erf}\left(\alpha-\beta\right)-1\right)$$
$$\frac{\partial f}{\partial\alpha}-2\alpha f(\alpha ,\beta)=\frac{\sqrt{\pi}}{2}\left(\exp\left(2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\exp\left(-2\alpha\beta\right)\ \operatorname{erf}\left(\alpha-\beta\right)-1\right)$$
$$\mu_2(\alpha)=\exp\left(-\int 2\alpha d\alpha\right)=\exp\left(-\alpha^2\right)$$
$$\left(\exp\left(-\alpha^2\right)f(\alpha,\beta)\right)'=\frac{\sqrt{\pi}}{2}\int\left(\exp\left(-\alpha^2+2\alpha\beta\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\exp\left(-\alpha^2-2\alpha\beta\right)\ \operatorname{erf}\left(\alpha-\beta\right)-\exp\left(-\alpha^2\right)\right)d\alpha$$
$$=-\frac{\operatorname{erf}\left(\alpha\right)\sqrt{\pi}}{2}+\frac{e\sqrt{\pi}}{2}\int\left(\underbrace{\exp\left(-\left(\alpha-\beta\right)^2\right)\ \operatorname{erfc}\left(\alpha+\beta\right)}_{IBP}-\exp\left(-\left(\alpha+\beta\right)^2\right)\ \operatorname{erf}\left(\alpha-\beta\right)\right)d\alpha$$
$$=-\frac{\operatorname{erf}\left(\alpha\right)\sqrt{\pi}}{2}+\frac{e\sqrt{\pi}}{2}\left(-\frac{\sqrt{\pi}}{2} \operatorname{erf}\left(\beta-\alpha\right)\ \operatorname{erfc}\left(\alpha+\beta\right)-\int\left(\underbrace{\operatorname{erf}\left(\beta-\alpha\right)+\operatorname{erf}\left(\beta+\alpha\right)}_{0}\right)\exp\left(-(\alpha+\beta)^2\right)d\alpha\right)$$
$$f(\alpha ,\beta)=-\frac{\sqrt{\pi}\operatorname{erf}\left(\alpha\right)e^{\alpha^2}}{2}-\frac{\pi e^{a^2+1}}{4}\left(\operatorname{erf}\left(\beta-\alpha\right)\ \operatorname{erfc}\left(\alpha+\beta\right)\right)+k$$
$$f(0,0)=\int_0^1\frac{dy}{1+y^2}=\frac{\pi}{4}=k$$
$$\boxed{f(\alpha,\beta)=-\frac{\sqrt{\pi}e^{\alpha^2}}{2}\operatorname{erf}\left(\alpha\right)-\frac{\pi e^{\alpha^2+1}}{4}\left(\operatorname{erf}\left(\beta-\alpha\right)\ \operatorname{erfc}\left(\alpha+\beta\right)\right)+\frac{\pi}{4}}$$
Therefore, the original integral is equal to: $$I=\exp\left(-\frac{a^2+b^2}{8T}\right)f\left(\frac{a}{2\sqrt{2T}},\frac{b}{2\sqrt{2T}}\right)=$$
$$ \bbox[5px,border:2px solid red] { \frac{\pi}{4}\exp\left(-\frac{a^2+b^2}{8T}\right)-\frac{\sqrt{\pi}}{2} \exp\left(-\frac{b^2}{8T}\right)\ \operatorname{erf}\left(\frac{a}{2\sqrt{2T}}\right)-\frac{\pi}{2} \exp\left(-\frac{b^2}{8T}+1\right)\ \operatorname{erf}\left(\frac{b-a}{2\sqrt{2T}}\right)\ \operatorname{erfc}\left(\frac{b+a}{2\sqrt{2T}}\right) } $$
Best Answer
Using the same approach
Expanded around $x=0$
$$\frac{1}{\sqrt{\tau }}\,\exp \Big[-\frac 12 (x-\tau )^2\Big] =\frac{1}{\sqrt{t}}\,\,e^{-\frac{t^2}{2}}\sum_{n=0}^\infty \frac {a_n}{n!}\,x^n$$ So,as you did, integrating with respect to $\tau$ from $0$ to $\infty$ to face $$v(x)=\sum_{n=0}^\infty b_n\, x^n$$ where $$b_n=\frac {2^{\frac{2n-3}{4} } } {n! }\,\,\Gamma \left(\frac{2n+1}{4} \right)\,\,\left(\sin \left(\frac{\pi n}{2}\right)+\cos \left(\frac{\pi n}{2}\right)\right)$$ is completely defined (as well as its derivatives.
This would be a fast convergent series since $$\left|\frac{b_{n+1}}{b_n}\right|\sim \frac 1{\sqrt n} \left(1-\frac{1}{n}\right)$$
Now, looking for the zero of $$v'_{(p)}(x)=\sum_{n=1}^p n\, b_n\, x^{n-1}$$
$$\left( \begin{array}{cc} p & x_{p} \\ 5 & 0.7957738274686 \\ 6 & 0.7693337827539 \\ 7 & 0.7627696942087 \\ 8 & 0.7645974040264 \\ 9 & 0.7651034151271 \\ 10 & 0.7649737593594 \\ 20 & 0.7649508673876 \\ \end{array} \right)$$
You should notice that, in principle, we do not need to solve the equation since, using power series reversion, we know its formal expression.