Suppose one is trying to find the first few terms of the asymptotic solution to an equation, which may be a differential equation, or a system of differential equations, in terms of a small parameter, $\delta$, in the limit $\delta\to0$.
Suppose that the equation has terms in $\delta^{-y}=e^{yln(1/\delta)}$ as well as terms in $\delta$, where $y$ is the dependent variable. One can imagine doing a series solution in powers of $\epsilon=1/ln(1/\delta)$: $$ y(x)=\sum_{n=0}^{\infty}y_{n}(x)\epsilon^{n} $$ Then, since $\epsilon^{n}>>\delta$, for all $n>0$, the $y_{1}(x)$ solution should be capturing the 'first order' behaviour. My question is, since this doesn't contain any corrections to deal with the $O(\delta)$ terms, how does one go about rectifying this? Should one do a double series expansion like: $$ y(x)=\sum_{n,m=0}^{\infty}y_{n}(x)\epsilon^{n}\delta^{m} $$ or is there a better way e.g. using a WKB approximation. I reckon there may be a general theory to describe this sort of problem but I haven't come across it. If anyone can point me in the right direction, that would be greatly appreciated.