I need to solve the second order equation
$$u^{\prime \prime}(t)=\frac{16 t\left(\beta t^3-27\right)}{81 \beta} u(t)$$
or alternatively,
$$u''(t) = \frac{16}{81}t^4u(t) - \frac{432}{81 \beta }tu(t)$$
$\beta > 0$. I'm not sure how to approach this as it is second order. The equation is in normal/regular form. Wolfram gives me a series solution, but I would like something closed form.
The original equation was
$$9f''-10t^2f'+\left(t^4+\left(\frac{48}{\beta}-10\right)t\right)f =0 $$
and I used the integrating factor
$$f(t) = \exp\left({\frac{5}{27}t^3}\right)u(t)$$
to obtain the form seen above.
EDIT: Here is a simpler related case which is solved:
$$9 f^{\prime \prime}-10 t^2 f^{\prime}+\left(t^4-10 t\right) f=0$$
which tranforms into regular form by $$ f(t)=\exp \left(\frac{5}{27} t^3\right) u(t) $$ where $u(t)$ then satisfies the equation $$ u^{\prime \prime}-\frac{16}{81} t^4 u=0 $$
the solutions of which are $$ t^{1 / 2} I_{ \pm 1 / 6}(T) \text { and } t^{1 / 2} K_{1 / 6}(T) $$ where $T=\frac{4}{27} t^3$ and $I_{ \pm 1 / 6}(T)$ and $K_{1 / 6}(T)$ are the modified Bessel functions of order $\frac{1}{6}$. In this way we see that the solutions of Eq. (1.4) are given by $$ t^{1 / 2} \exp \left(\frac{5}{27} t^3\right)\left\{I_{ \pm 1 / 6}(T), K_{1 / 6}(T)\right\}. $$
In the original equation, $$9 f''(t) - 10 t^2 f'(t) + \left(t^4 + \left(\frac{48}{\beta} - 10\right)\right) f(t) = 0 ,$$ changing variables via $f(t) = t \exp\left(\frac{t^3}{27}\right) w(t)$ and $\tau = \frac{8 t^3}{27}$ gives $$\tau w''(\tau) + \left(\frac{4}{3} - \tau\right) w'(\tau) - \left(\frac{2}{3} - \frac{2}{\beta}\right) w(\tau) = 0,$$ which is Kummer's equation with parameters $\mu = \frac{2}{3} - \frac{2}{\beta}, \nu = \frac43$, so the general solution is $$w(\tau) = c_1 M_{\frac{2}{3} - \frac{2}{\beta}, \frac43}(\tau) + c_2 U_{\frac{2}{3} - \frac{2}{\beta}, \frac43}(\tau) ,$$ where $M_{\mu, \nu} = {}_1 F_1(\mu; \nu; \,\cdot\,)$ is Kummer's confluent hypergeometric function, and $U_{\mu, \nu}$ is Tricomi's confluent hypergeometric function; for more, see the above link, as well as the NIST DLMF entry for Kummer functions.
Translating back to $f(t)$ thus gives the general solution $$\boxed{f(t) = c_1 t \exp \left(\frac{t^3}{27}\right) M_{\frac{2}{3} - \frac{2}{\beta}, \frac43}\left(\frac{8t^3}{27}\right) + c_2 t \exp \left(\frac{t^3}{27}\right) U_{\frac{2}{3} - \frac{2}{\beta}, \frac43}\left(\frac{8t^3}{27}\right)} .$$
Remark This form of the solution is useful in part because it lets us readily identify parameter values $\beta$ and limiting cases that yield special behavior.
When $\beta = -3$, we have $\mu = \frac23 - \frac2\beta = \frac43 = \nu$, but $M_{\mu, \mu} = \exp$ and $U_{\mu, \mu}(\tau) = e^\tau \Gamma(1 - \mu, \tau)$, where $\Gamma$ is the incomplete gamma function, so the general solution in $w$ is $$w \tau = c_1 e^\tau + c_2 \Gamma\left(-\frac13, \tau\right) .$$
When $\beta = 3$, the equation in $w$ reduces to first order in $y := w'$, $$\tau y'(\tau) + \left(\frac43 - \tau\right) y(\tau) = 0 ,$$ which has general solution $$y(\tau) = C \tau^{-\frac43} e^\tau ,$$ so in this case the solutions $w(\tau)$ are $$w(\tau) = c_1 + c_2 \int^\tau \upsilon^{-\frac43} e^\upsilon \,d\upsilon ,$$ the second terms of which can be expressed in terms of the incomplete gamma function.
When $\beta = 6$, we have $\mu + 1 = \nu$, and $U_{\mu, \mu + 1}(\tau) = \tau^{-\mu}$; on our case $\mu = -\frac23$, which gives the general solution $$w(\tau) = c_1 \tau^{-\frac13} + c_2 \int^\tau \upsilon^{-\frac23} e^\upsilon \,d\upsilon .$$
When $\beta = \frac{6}{2 + 3 k}$, $k \in \{1, 2, 3, \ldots\}$, $\mu = -k$, and $M_{-k, \nu}$ is a generalized Laguerre polynomial of degree $k$ (and $M_{-k, \mu}$ and $U_{-k, \mu}$ are linearly dependent). For example, for $k = 1, 2$ ($\beta = \frac65, \frac34$, respectively), the general solutions are $$ \begin{array}{ccr} \hline k & \beta \qquad & w(\tau) \qquad \qquad \qquad \qquad \\ \hline 1 & \frac65 \qquad & c_1\left(1 - \frac34 \tau\right) + c_2 \tau^{1 / 3} M_{-\frac43, \frac23}(\tau) \\ 2 & \frac34 \qquad & c_1\left(1 - \frac32 \tau + \frac9{28} \tau^2\right) + c_2 \tau^{1 / 3} M_{-\frac73, \frac23}(\tau) \\ \hline \end{array} .$$
Finally, as you've already pointed out, in the limit $\beta \to \pm\infty$, we can express the solution in terms of (modified) Bessel functions: The parameter values are related by $2 \mu = \nu$, and in this case $M_{\mu, \nu}$ satisfies the identity $$M_{\mu, 2 \mu}(\tau) = 2^{2 \mu - 1} \tau^{\frac12 - \mu} e^{\frac\tau2} \Gamma\left(\mu + \frac12\right) I_{\mu - \frac12} \left(\frac\tau2\right) ,$$ and in our case ($\mu = \frac23$), $$M_{\frac23, \frac43}(\tau) = \sqrt[3]{2} \tau^{-\frac16} e^{\frac\tau2} \Gamma\left(\frac76\right) I_{\frac16} \left(\frac\tau2\right) ,$$ recovering the appearance of the Bessel function $I$ of order $\frac16$. There's a similar identity relating $U_{\mu, 2 \mu}$ and the Bessel function $K_{\mu - \frac12}$.