Laplace's Method to Arbitrarily Higher Order

This implements Laplace's method for the integral  (∫_a^b)(e^((1/g)S[ϕ]))dϕ  when S has a unique maximum  at some point  a<ϕmax<b.

We can expand
                   S[ϕ]=S[ϕmax]-(1/(2G))+(∑_(k=3))(Subscript[S, k])(((ϕ-ϕmax))^k)
with G>0. Then 
                              (∫_a^b)(e^((1/g)S[ϕ]))dϕ =(e^((1/g)(Subscript[S, 0])))(Sqrt[2ℼg G])Y[g]
and
                                          Y[g]~1+(∑_r)(g^r)(Subscript[Y, r]). 
This program  this  asymptotic expansion to any given  order in g,  limited by the memory of the computer. 

We verify the answer to the 11th order against the built-in asymptotics of two standard examples: the factorial and the BesselK[1] function.

How do you tell Mathematica that (0^(0 ))=1? Currently using an awkward workaround by using a function power[x,n].

Clear[ϕ,S,ϕmax,G,s,power] (*Clear some quantities that will be used later*)

Term[order_]:=Module[{kmax},
kmax=2*(order+1);
weights=Table[(1/2)k-1,{k,3,kmax}];
parts=IntegerPartitions[order,2*order,weights];
(*Find all the vertices and their multiplicities at this order*)
numterms=Dimensions[parts][[1]];
ks=Table[2*parts[[i]]+2,{i,1,numterms}];
ns=Table[Table[Count[ks[[i]],k],{k,3,kmax}],{i,1,numterms}];
edges=Table[(1/2)Sum[(j+2)*ns[[i,j]],{j,1,kmax-2}],{i,1,numterms}];


Sum[(2*edges[[i]]-1)!!(G^(edges[[i]]))
,{i,1,numterms}];
(*Table[((s[k])^(ns[[1,kmax-2]])),{k,3,kmax}]*)
(g^order)*Sum[(2*edges[[i]]-1)!! power[G,edges[[i]]]
Product[((power[s[k],ns[[i,k-2]]])/(ns[[i,k-2]]!)),{k,3,kmax}],{i,1,numterms}]
]
Laplace[ordermax_]:=Sum[Term[order],{order,0,ordermax}]

Explicit answer to order 4

power[x_,n_]=x^n;

Laplace[4]

Verifying the formula to the 11th order  in the Stirling  approximation

Clear[ϕ,S,ϕmax,G,s,power]
power[x_,n_]=
If[(n==0&&x==0),1,(x^n)];
S[ϕ_]=-ϕ+Log[ϕ]
ϕmax=1;
S[ϕmax];
s[k_]:=Limit[((1/(k!)))D[S[ϕ],{ϕ,k}],ϕ->ϕmax]
(*s[k_,ϕ_]=((1/(k!)))D[S[ϕ],{ϕ,k}];*)

G=((-1)/(2s[2]));
Laplace[11]
Series[((n!)/((n^(n+(1/2)))Exp[-n](Sqrt[2Pi]))),{n,Infinity,11}]

Verifying with the BesselK Function of order one

Clear[ϕ,S,ϕmax,G,s,power]
power[x_,n_]=
If[(n==0&&x==0),1,(x^n)];
S[ϕ_]=((-1)/2)(ϕ+(1/ϕ));
ϕmax=1;


s[k_]:=Limit[((1/(k!)))D[S[ϕ],{ϕ,k}],ϕ->ϕmax];
s[4];

G=((-1)/(2s[2]));
Laplace[11]
Refine[Integrate[Exp[((-1)/(2g))(t+(1/t))],{t,0,Infinity}],Assumptions->(g>0)]

Refine[Series[(((Exp[(1/g)S[ϕmax]](Sqrt[2Pi*G*g])))^(-1))2BesselK[1,(1/g)],{g,0,11}],Assumptions->(g>0)]