2nd order linear ODE non constant coefficient solution methods I'm working on a project and I've come across a seemingly standard ODE, however I have no idea how to solve it.  The equation in question is 
$$
\epsilon'' + \left[g(\tau) + \frac{\beta}{\tau}\right]\epsilon' + \left[\frac{\gamma}{\tau^{4/3}}\right]\epsilon = 0
$$
where $g(\tau) = \alpha \cos(\tau)$ or $\delta\cos(\tau)\sin(\tau)$ (there are two cases).  I've tried a series solution but unsurprisingly it's pretty messy and doesn't really solve anything, and I've tried Laplace transforms but unfortunately they're divergent for negative powers of $\tau$.  
Does anybody have recommendations for methods I can try, or even how to make a decent ansatz?  This is just a step in my project and it's really presenting a roadblock for me, so any help is appreciated.  
 A: You have a linear time-variant differential equation. You can cast it into matrix form by introducing $\epsilon = \epsilon_1$, $\epsilon'=\epsilon_1'=\epsilon_2$. You will get
$$\dfrac{d}{d\tau}\begin{bmatrix} \epsilon_1\\\epsilon_2 \end{bmatrix}=\begin{bmatrix} 0 & 1\\-\dfrac{\gamma}{\tau^{4/3}} & -g(\tau)-\dfrac{\beta}{\tau} \end{bmatrix}\begin{bmatrix} \epsilon_1\\\epsilon_2 \end{bmatrix}.$$
In this form, you can invoke the Peano-Baker series to obtain the solution
$$\boldsymbol{\epsilon}(\tau)=\left[\boldsymbol{I}+\int_{\tau_0}^{\tau}\boldsymbol{A}(\tau_1)d\tau_1+\int_{\tau_0}^\tau\int_{\tau_0}^{\tau_1}\boldsymbol{A}(\tau_1)\boldsymbol{A}(\tau_2)d\tau_2d\tau_1+... \right]\boldsymbol{\epsilon}(\tau=\tau_0),$$
in which 
$$\boldsymbol{A}(\tau)=\begin{bmatrix} 0 & 1\\-\dfrac{\gamma}{\tau^{4/3}} & -g(\tau)-\dfrac{\beta}{\tau} \end{bmatrix},$$
$\boldsymbol{I}$ is the $2\times 2$ identity matrix, and $\boldsymbol{\epsilon}=[\epsilon_1,\epsilon_2]^T$.
If there is a closed form solution you should be able to retrieve it with this method.
