I have the following recurrence
$$
\left(4-3 n^3\right) \mu _{n+1}+n ((n-3) n+3) \mu _n+(3 n (n (n+3)+3)+4) \mu _{n+2}+6\\
=(n+2)^3 \mu _{n+3}
$$
with initial values
$$
\left\{0,1,\frac{9}{8},\frac{281}{216},\frac{9955}{6912},\frac{1340387}{864000},\frac{8519417}{5184000},\frac{3061810421}{1778112000}\right\}$$
starting from $n=1$.
I have tried to use Sigam to solve it. But it can only find one inhomogenious solution
$$
\frac{6}{11}\left(\underset{\iota _1=2}{\overset{n}{\sum}}\frac{1}{\iota _1-1}\right)
$$
Is it possible to find an homogenious solution?
If not possible, how can we get an asymptotic expansion of $\mu_n$?
I am able to estimate numerically that
$$
\mu_n \approx \frac{6}{11}\left(\underset{\iota _1=2}{\overset{n}{\sum}}\frac{1}{\iota _1-1}\right) + 0.308492+0.0565532 \left(n^{-3+i \sqrt{2}}+n^{-3-i \sqrt{2}}\right)
$$
using the Asymptotics package. I wonder if it is possible to find a closed form for constant term.
Best Answer
The constants you sought may be expressed in terms of (generalized) hypergeometric functions. As the results are quite involved, I shall only sketch the procedure.
Put $$\mu_n=\nu_n+\frac6{11}\sum_{\iota=1}^{n-1}\frac1\iota,$$ $\nu_n$ satisfies the homogeneous recurrence. Let $$f(z)=\sum_{n=1}^\infty \nu_nz^n$$ be the generating function. The recurrence translates to the inhomogeneous DE $$z^3 (1-z)^3 f'''(z)+z (1-z)^3 f'(z)-(1-2 z+7 z^2) f(z)=z^2 (u+4vz),$$ where $$\begin{aligned} &u=\nu_2-\nu_1,&&v=2\nu_3-\nu_2-\nu_1 \end{aligned}$$ depend on initial values. The singularities of this DE lie at $z=0$, $1$, $\infty$. The asymptotic behavior of $\nu_n$ is reflected on the singularity of $f(z)$ at $z=1$. At this point the DE has a holomorphic particular integral $P(z)$, and a set of fundamental solutions of the homogeneous DE are $$\begin{aligned} &\frac{z}{1-z},&&F_1(z)=(1-z)^{2+i\sqrt{2}}(1+c_1(1-z)+\cdots),&&F_2(z)=(1-z)^{2-i\sqrt{2}}(1+c_2(1-z)+\cdots). \end{aligned}$$ If we have $$f(z)=P(z)+C\frac{z}{1-z}+C_1F_1(z)+C_2F_2(z),$$ then the standard procedure (Darboux's method) yields $$\nu_n\sim C+C_1\frac{\sinh (\sqrt{2} \pi ) \Gamma (3+i \sqrt{2})}{\pi i}n^{-3-i\sqrt{2}}-C_2\frac{\sinh (\sqrt{2} \pi ) \Gamma (3-i \sqrt{2})}{\pi i}n^{-3+i\sqrt{2}}.$$ In general the coefficients $C$'s can only be determined numerically, as they depend on the special values of these solutions.
In the present case, however, a reduction to the hypergeometric DE is possible. Put $$f(z)=\frac{z}{1-z}\int g(z)dz,$$ it leads to $$z^2 (1-z)^2g''(z)+3 z (1-z) g'(z)+(1+4 z+z^2)g(z)=u+4vz.$$ If we further substitute $$g(z)=z^{-1}(1-z)^{2+i\sqrt{2}}h(z),$$ we find $$z (1-z) h''(z)+(1-2(1+i \sqrt{2}) z) h'(z)+(1-i \sqrt{2}) h(z)=(1-z)^{-3-i\sqrt{2}}(u+4vz).$$ The left hand side is hypergeometric DE, so we can find a set of fundamental solutions of the homogeneous DE of $g(z)$ $$G_1(z)=\frac{(1-z)^{2+i \sqrt{2}}}{z}\, _2F_1\left(\frac{1}{2}+i \sqrt{2}-\frac{i \sqrt{3}}{2},\frac{1}{2}+i \sqrt{2}+\frac{i \sqrt{3}}{2};1+2 i \sqrt{2};1-z\right),$$ $$G_2(z)=\frac{(1-z)^{2-i \sqrt{2}}}{z}\, _2F_1\left(\frac{1}{2}-i \sqrt{2}-\frac{i \sqrt{3}}{2},\frac{1}{2}-i \sqrt{2}+\frac{i \sqrt{3}}{2};1-2 i \sqrt{2};1-z\right).$$ By variation of constants, we find $$\begin{aligned} g(z)&=G_1(z)\int_0^z\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_2(z)dz-G_2(z)\int_0^z\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_1(z)dz\\ &=Q(z)+c_1G_1(z)+c_2G_2(z), \end{aligned}$$ where $$Q(z)=-G_1(z)\int_z^1\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_2(z)dz+G_2(z)\int_z^1\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_1(z)dz$$ is holomorphic at $z=1$, and $$\begin{aligned} &c_1=\int_0^1\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_2(z)dz,&&c_2=-\int_0^1\frac{iz(u+4vz)}{2\sqrt{2}(1-z)^5}G_1(z)dz. \end{aligned}$$ Note that these integrals actually diverges, but this difficulty can be circumvented by analytic continuation or by other means.
Finally, $$f(z)=\frac{z}{1-z}\int_0^zg(z)dz=\frac{z}{1-z}\int_0^1g(z)dz-\frac{z}{1-z}\int_z^1(Q(z)+c_1G_1(z)+c_2G_2(z))dz.$$ It follows that $$\begin{aligned} &C=\int_0^1g(z)dz,&&P(z)=-\frac{z}{1-z}\int_z^1Q(z)dz,&&C_iF_i(z)=-c_i\frac{z}{1-z}\int_z^1G_i(z)dz. \end{aligned}$$ In particular, $$\begin{aligned} &C_1=\frac{-3+i \sqrt{2}}{11}c_1,&&C_2=\frac{-3-i \sqrt{2}}{11}c_2. \end{aligned}$$
In summary, the second constant 0.0565532 can be expressed by integrals involving $\, _2F_1$, or by values of $\, _3F_2$ at $1$. The first constant 0.308492 is expressed by double integrals involving $\, _2F_1$, which I do not know whether it can be reduced further or not.
I retain the dependence on initial values for there is a special case that renders a more exact treatment. If $u=4v=3k$ or $\nu_2=\nu_1+3k$, $\nu_3=\nu_1+15k/8$, the DE of $g(z)$ has a particular solution $Q(z)=k/z$, and the constants in the asymptotic expansion can be expressed by Gamma functions.