Asymptotic expansion of $\sum\limits_{n=1}^{\infty} \frac{x^{2n+1}}{n!{\sqrt{n}} }$

I sketch the arguments for $C(x)$, the arguments for $L(x)$ are essentially the same.

The specific form of the sum suggests probabilistic arguments. Let $X_x$ be a $\mathrm{Poiss}(x^2)$-distributed random variable and note that $$C(x)\,e^{-x^2}=\mathbb{E}\frac{x}{\sqrt{X_x}}\,1_{\{X_x\geq 1\}}$$ It is known that $\frac{X_x -x^2}{x}\longrightarrow N(0,1)$ (standard normal) in distribution as $x\longrightarrow \infty$ and that that all moments $\mathbb{E}\left(\frac{X_x-x^2}{x}\right)^k$ of $X_x$ converge to the corresponding moments of $N(0,1)$, so that (for large $x$) $X_x$ is concentrated around $x^2$, with deviations of order $x$.

To use that information split $\{X_x\geq 1\}$ on the rhs into (say) the parts $1\leq X_x <\tfrac{1}{2}x^2$, $X_x-x^2 >\tfrac{1}{2} x^{2}$ and $|X_x-x^2| \le \tfrac{1}{2}x^{2}$.

By routine arguments the integrals over the first two parts are asymptotically exponentially small. For the remaining part write \begin{align*} \mathbb{E}\frac{x}{\sqrt{X_x}}\,1_{\{|X_x-x^2|\leq \tfrac{1}{2}x^2\}} &=\mathbb{E}\frac{x}{\sqrt{x^2+(X_x-x^2)}}\,1_{\{|X_x-x^2|\leq \tfrac{1}{2} x^2\}}\\ &=\mathbb{E}\frac{1}{\sqrt{1 +\frac{1}{x^2}(X_x-x^2)}}\,1_{\{|X_x-x^2|\leq \tfrac{1}{2} x^2\}}\\ &=\mathbb{E}\sum_{k=0}^\infty {-\frac{1}{2} \choose k} \frac{1}{x^{2k}}(X_x-x^2)^k \,1_{\{|X_x-x^2|\leq \tfrac{1}{2}x^2\}}\\ %C(x)\,e^{-x^2}&= 1+\frac{3}{8}x^{-2} + \frac{65}{128}x^{-4} + \frac{1225}{1024}x^{-6} + \frac{1619583}{425984}x^{-8}+\mathcal{O}(x^{-10})\\ %L(x)\,x^2\,e^{-x^2}&= 1+\frac{15}{8}x^{-2} + \frac{665}{128}x^{-4}+\frac{19845}{1024}x^{-6}+\frac{37475823}{425984}x^{-8}+\mathcal{O}(x^{-10}) \end{align*} Clearly the series may be integrated termwise, and completing the tails changes it only by asymptotically exponentially small terms. For $k\geq 2$ the central moment $c_k(x):=\mathbb{E}\left(X_x-x^2\right)^k$ is a polynomial in $x^2$ of degree $\lfloor k/2\rfloor$. Thus (after regrouping of terms) the formal series $$\sum_{k=0}^\infty x^{-2k}{-\tfrac{1}{2} \choose k} c_k(x)$$ gives a full asymptotic expansion of $C(x)e^{-x^2}$. Evaluating the first eight terms gives $$C(x)\,e^{-x^2}= 1+\frac{3}{8}x^{-2} + \frac{65}{128}x^{-4} + \frac{1225}{1024}x^{-6} + \frac{131691}{32768}x^{-8}+\mathcal{O}(x^{-10})$$

EDIT: I corrected the coefficient of $x^{-8}$. Thanks to Johannes Trost for pointing out that I miscalculated. Similarly the coefficient of $x^{-8}$ in the asymptotic series for $L(x)$ given in the comment must be corrected.

Warning: Nearly every number in this answer is wrong ! Please, read the answer posted by esg, which gives the right asymptotic expansion !

Such problems can be solved by Laplace's method. The starting point is the observation, that the values of the terms of the sums are unimodal as function of $n$ for fixed (large) $x$, i.e., have one maximum. If the maximum (as in the given cases) is sharp enough, the terms of the sum can be approximated by a Gaussian. This Gaussian is then integrated on the whole $n$-axis, where $n$ is considered as a real number. The integral often converges. The contributions from negative $n$ are negligible.

The maximum, $n_{0}$, for $\frac{x^{2 n +1}}{n!\sqrt{n}}$ for fixed $x$ can be evaluated to $$ n_{0} = x^2-1-\frac{5}{12}x^{-2}-\frac{1}{2} x^{-4}-\frac{123}{160} x^{-6}+\frac{359}{180} x^{-8} + O(x^{-10}). $$ The sum is approximated by an integral over a Gaussian (which we get when expanding around $n=n_{0}$ to second order in $\tau$) $$ C(x) \approx \int_{-\infty}^{\infty} d\tau\ \exp\left( \ln\left(\frac{x^{2 n +1}}{n!\sqrt{n}}\right)|_{n\rightarrow n_{0}+\tau}\right), $$ consistently expanded to $O(x^{-10})$ after integration. The result is $$ C(x)= e^{x^2}\left(1+\frac{5}{12} x^{-2}+\frac{157}{288} x^{-4}+\frac{49729}{51840} x^{-6}+\frac{417001}{497664} x^{-8}+O(x^{-10})\right). $$ I could not find any systematic for the coefficients. Only the numerators of the higher order terms seem to contain rather large prime factors.

The result is different from your findings, though. However, numerical evidence suggests that the above asymptotic expansion is correct.

My result for $L(x)$ is $$ L(x)=e^{x^2}\left(x^{-2}+\frac{23}{12} x^{-4}+\frac{1525}{288} x^{-6}+\frac{949099}{51840} x^{-8}+O(x^{-10})\right), $$ with again rather large prime factors of the nominators of the higher order coefficients. The maximum of the terms as function of $n$ is reached (up to $O(x^{-10})$) for $n=n_{0}$ with $$ n_{0}=x^2 -2-\frac{23}{12}x^{-2}-5 x^{-4}-\frac{2643}{160}x^{-6}-\frac{3007}{45} x^{-8}+O(x^{-10}). $$ A numerical test shows very good quality of the asymptotic expansions for $C(x)$ and $L(x)$, even for $x$ around $2$.

All calculations where done with Mathematica 11.

Edit: I corrected a typo in the coefficient of $x^{-8}$ in the expansion of $n_{0}$ for the asymptotics of $C(x)$.

Edit: More numerical calculations indicate that the asymptotic expansion for $C(x)$ given by me is by far not as accurate as I hoped. The numerical evidence shows that $$ C(x) = e^{x^{2}}\left(1+\frac{3}{8}x^{-2}+\frac{1}{2}x^{-4}...\right) $$ is much better.

Edit: The coefficients for $L(x)$ as given in my answer are wrong as well. The numbers given by esg in the comment below seem to be correct.

This is another, totally different (and correct !) approach for answering the question. It is simply too long for a comment. So I decided to write it in a new answer. (Although that might look odd, but my previous answer, though accepted, is wrong.)

First define $$ I_{\mu}(y) = y^{1/2} \sum_{n=1}^{\infty} \ \frac{y^{n}}{n! \ n^{\mu}}, $$ and observe that the OP's functions are $$ C(x)= I_{\frac{1}{2}}(x^{2}) $$ and $$ L(x)=I_{\frac{3}{2}}(x^{2}) $$ Let $m$ be the integer part of $\mu$ and $\xi$ the fractional part, $0<\xi<1$.

Replace $n^{-\mu}$ in the definition of $I_{\mu}(y)$ by a ratio of $\Gamma$-functions times an asymptotic series in $n$ using this formula and for the coefficients of the asymptotic series in $n$ using the Norlund polynomials, $B_k^{(\alpha)}(x)$ found here. They are available in Mathematica via $\mathtt{NorlundB[k,\alpha,x]}$. Concretely, $$ n^{-\xi}\sim \frac{\Gamma(n)}{\Gamma(n+\xi)}\sum_{k=0}^{\infty}n^{-k} {\xi \choose k} B_k^{(1+\xi)}(\xi). $$

Now exchange the summation of the asymptotics in $n$ (with, say, summation index $k$) and the summation over $n$, which I assume light heartedly to be possible. The (now inner) sum over $n$ results in generalized hypergeometric functions of the form $$ _{m+k+2}F_{m+k+2}\left(\left. {1,\dots,1}\atop{2,\dots,2,1+\xi} \right\vert y\right), $$ with $m+k+2$ 1s in the upper line and $k+m+1$ 2s in the lower line. To get there one has to shift the summation index such that summation starts at $n=0$. Then insert $n+1=\frac{(2)_{n}}{(1)_{n}}$, with the usual Pochhammer symbols used. The defining formula for generalized hypergeometric functions results.

Using the asymptotic expansion of the generalized hypergeometric function for $y\rightarrow \infty$ given in a paper by Volkmer and Wood (downloadable from here) and after some simplifications one arrives at the asymptotic formula $$ I_{m+\xi}(y)=e^{y}\ y^{\frac{1}{2}-m-\xi}\ \frac{\xi\ \sin \pi\xi}{\pi}\sum_{k=0}^{\infty}(-1)^{k+1} y^{-k} \frac{\Gamma(k-\xi)}{k!}\ B_{k}^{(1+\xi)}(\xi) \\ \sum_{s=0}^{\infty} y^{-s} \left\{ {\sum_{(s_{1},\dots,s_{m+k+1})}} \frac{\Gamma(\xi + s_{m+k+1})}{s_{m+k+1}!}\prod_{j=1}^{m+k+1} \frac{\Gamma\left(j+\sum_{i=1}^{j}s_{i}\right)}{\Gamma\left(j+\sum_{i=1}^{j-1}s_{i}\right)}\right\}. $$ $(s_{1},\dots ,s_{m+k+1})$ under the sum sign indicates summation over all (ordered) partitions of $s$ into $m+k+1$ non-negative integers, $s_{1},\dots ,s_{m+k+1}$. Order matters here, i.e., $(1,0)$ is different from $(0,1)$.

Numerical calculations of the coefficients (with highest possible precision on my laptop) show excellent match (5 or more digits) with the formula, even for exotic indices, like $\mu=\pi$ and orders up to $y^{-6}$.

For the OP's functions I get: $$ C(x) e^{-x^{2}} = 1 + \frac{3}{8} x^{-2} + \frac{65}{128} x^{-4} + \frac{1225}{1024} x^{-6} + \frac{131691}{32768} x^{-8} + O(x^{-10}) , $$ $$ L(x) x^{2} e^{-x^{2}} = 1 + \frac{15}{8} x^{-2} + \frac{665}{128} x^{-4} + \frac{19845}{1024} x^{-6} + \frac{2989371}{32768} x^{-8} + O(x^{-10}) . $$

All calculations were done with Mathematica 11.