Boundary layer in time Consider the initial value problem $\varepsilon x'' + x' + tx = 0$ where $x(0) = 0$ and $x'(0) = 1$. I'm solving this problem using a matched asymptotic expansion. First, I let $$x(t, \varepsilon) = \varepsilon x_1(t) + \varepsilon^2 x_2(t) + o(\varepsilon^2),$$ and solve at corresponding orders of $\varepsilon$ and leave the coefficients undetermined. Then I define $T = t/\varepsilon$ and let $X(\varepsilon, T) = x(\varepsilon, t)$ and write the re-scaled equation $$X'' + X' + \varepsilon TX = 0.$$
Then solve this at the corresponding orders of $\varepsilon$ for the expansion $$X(\varepsilon, T) = \varepsilon X_1(T) + \varepsilon^2 X_2(T) + o(\varepsilon^3).$$ Assuming $X$ is my solution 'inside' the boundary layer, I let $X$ satisfy the initial conditions $X(0) = 0, X'(0) = \varepsilon$. Then $X_1(0) = 0, X_1'(0) = 1$ and $X_n(0) = 0, X_n'(0) = 0$ where $n>1$. However, when I solve for the solutions 'inside' the boundary layer I get imaginary error functions. I believe this will cause me problems later in the matching.
I believe I have solved this problem without the $t$ term and I am building up to a more complicated function. I am familiar with boundary laters in the context of space, but cannot find much on these 'initial layer' problems.
I know I have not provided my solutions here, but I want to know if my method is sensible or if someone knows a better method of solving this (equations of this type). I appreciate your response.
 A: There appears to be a scale error in your re-balanced equation that might throw off later computations.
Perturbation series
Note: As $x''(0)=-ε^{-1}$, etc., the first terms of the Taylor expansion are $x-\frac12ε^{-1}x^2+O(x( ε^{-1}x)^2)$ which has a peak at $x=ε$ of magnitude $ε/2$. This means that simultaneously to rescaling the time, it makes sense to also scale to compensate the function so that the initial slope remains $1$.
With $X(T)=ε^{-1}x(εT)$ and thus $X'(T)=x'(εT)$, $X''(T)=εx''(εT)$ you should get
$$
X''(T)+X'(T)=εx''(εT)+x'(εT)=-εTx(εT)
\\
\implies X''(T)+X'(T)+ε^2TX(T)=0, ~~X(0)=0,~~X'(0)=1.
$$
which means that the perturbation parameter is $ε^2$, $X(T)=X_0(T)+ε^2X_1(T)+ε^4X_2(T)+...$.
The first order approximation is $X_0(T)=1-e^{-T}$.  The next term is obtained via
\begin{align}
X_1''(T)+X_1'(T)&=-TX_0(T)=-T+Te^{-T}, ~~X_1(0)=0,~~X_1'(0)=0\\
X_1(T)&=-\tfrac12T^2-(\tfrac12T^2+T+1)e^{-T}
\end{align}
Further terms will have higher polynomial degrees, leading to divergence for $T\to \infty$ and thus no match to the outer solution.
However, knowing that the outer solution is $Ce^{-t^2/2}=Ce^{-ε^2T^2/2}=C(1-\frac12ε^2T^2+...)$, one can absorb the obvious terms of the perturbation expansion in similar terms
\begin{align}
X(T)=X_0(T)+ε^2X_1(T)+...&=1-\tfrac12T^2-(1+ε^2(\tfrac12T^2+T+1))e^{-T}+...\\
&=e^{-ε^2T^2/2}-e^{ε^2T^2/2}e^{-T} - ε^2(T+1)e^{-T}
\end{align}
or
$$
x(t)=εe^{-t^2/2}-εe^{t^2/2}e^{-t/ε}-ε^2(t+ε)e^{-t/ε}
$$
(or possibly also $x(t)=εe^{-t^2/2}-ε(2-e^{-t^2/2})e^{-t/ε}+...$, depending on the higher order terms)
But this is guess-work that could be invalidated with every new term in the perturbation series.

Two time-scale approach
Coming from a slightly different angle, starting from the inner solution, set the integration constants there directly as "slowly moving" functions of $t$, that is, try to find a two-scale solution as $$x(t)=εA(t)-εB(t)e^{-t/ε}.$$ Then immediately $A(0)=B(0)$. With the derivatives one finds the other initial condition and the insertion into the differential equation.
\begin{align}
x'(t)&=εA'(t)+(B(t)-εB'(t))e^{-t/ε}, \\
εx''(t)&=ε^2A''(t)-(B(t)-2εB'(t)+ε^2B''(t))e^{-t/ε}, \\ \hline
0=εx''(t)x'(t)+tx(t)&=ε[εA''(t)+A'(t)+tA(t)] - ε[-B'(t)+εB''(t)+tB(t)] e^{-t/ε},\\
1&=ε(A'(0)-B'(0))+B(0).
\end{align}
In first order, separating the terms in $A$ and $B$ and using only the lowest order terms in $ε$, one finds the coefficients as $A_0(t)=e^{-t^2/2}$ and $B_0(t)=e^{t^2/2}$. 
In the next order of $A(t)=A_0(t)+εA_1(t)$, $B(t)=B_0(t)+εB_1(t)$,
$$
(e^{t^2/2}A_1(t))'=-(t^2-1)\implies A_1(t)=(t-\frac13t^3)e^{-t^2/2}\\
(e^{-t^2/2}B_1(t))'=(t^2+1)\implies B_1(t)=(t+\frac13t^3)e^{t^2/2}
$$
so that $A'(0)=-ε$, $B'(0)=ε$, and again $A(0)=B(0)=1$
Plotting plots of these first and second order approximations against the numerical solution gives a good fit.

A: WKB approximation
Look for basis solutions of the form $x(t)=\exp(S(t)/ε)$. Then $εx'(t)=S'(t)\exp(S(t)/ε)$ and $ε^2x''(t)=[εS''(t)+S'(t)^2]\exp(S(t)/ε)$. Inserting and canceling the exponential gives
$$
0=e^{-S/ε}(ε^2x''+εx'+εtx)=εS''(t)+S'(t)^2+S'(t)+εt
\\~~\\
\iff S'(t)^2+S'(t)=-ε(S''(t)+t).
$$
For simplicity name $s(t)=S'(t)$ and compute the terms of the perturbation series $s(t)=s_0(t)+εs_1(t)+...$
\begin{array}{rlrl|rl}
s_0^2+s_0&=0                     &\implies s_0&=0      &\text{ or }~~ s_0&=-1\\
2s_0s_1+s_1&=-t                  &\implies s_1&=-t     & s_1&=t \\
s_1^2+2s_0s_2+s_2&=-s_1'         &\implies s_2&=1-t^2  & s_2&=1+t^2\\
\end{array}
The approximation so far is
$$
x(t)=A\exp(-\tfrac12t^2+ε(t-\tfrac13t^3))+B\exp(-ε^{-1}t+\tfrac12t^2+ε(t+\tfrac13t^3))
$$
with $0=x(0)=A+B$ and $1=x'(0)=-ε^{-1}B$, so that $A=ε$, $B=-ε$.
The plot of these two approximations against the numerical solution gives a good fit even for largish values of $ε$.

A: This is not an answer to the question, but I think it might be useful nevertheless.
This equation can be solved exactly in terms of Airy functions. Let $x(t)=e^{\lambda t}W(t)$ and substitute into the differential equation to give (after factoring out $e^{\lambda t}$,
$$\epsilon\lambda^2 W(t)+2\epsilon\lambda W'(t)+\epsilon W''(t)+\lambda W(t)+W'(t)+tW(t)=0. $$
If $2\epsilon\lambda+1=0$ then we can remove the $W'$ terms, so $\lambda=-1/(2\epsilon)$ and $u=e^{-t/(2\epsilon)}$. This represents the slow decay in the solution for $x(t)$.
Now we are left with $$\frac{1}{4\epsilon} W+\epsilon W''-\frac{1}{2\epsilon}W+tW=0\Rightarrow \epsilon^2 W''+\left(\epsilon t-\frac{1}{4}\right)W=0.$$
The initial conditions are, in terms of $W$, $W(0)=0$ and $W'(0)=1$. Now let $$s=\frac{1-4\epsilon t}{4\epsilon^{4/3}}$$ and $H(s)=W(t)$ so that
$ H'(s)=-W'(t)\epsilon^{-1/3}$ and $H''(s)=\epsilon^{-2/3}W''(t)$. Then
$$H''(s)-sH(s)=0,\quad H\left(\frac{1}{4\epsilon^{4/3}}\right)=0,\quad H'\left(\frac{1}{4\epsilon^{4/3}}\right)=-\epsilon^{1/3}.$$
This is the Airy differential equation, and its solution is a combination of Airy functions $\textrm{Ai}(s)$ and $\textrm{Bi}(s)$,
$$ H(s) =  c_1\textrm{Ai}(s)+c_2\textrm{Bi}(s),$$
and $c_1$ and $c_2$ satisfy
$$ c_1\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)+c_2\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)=0,\quad c_1\textrm{Ai}'\left(\frac{1}{4\epsilon^{4/3}}\right)+c_2\textrm{Ai}'\left(\frac{1}{4\epsilon^{4/3}}\right)=-\epsilon^{1/3}, $$
or 
$$c_1 = \epsilon^{1/3}\frac{-\textrm{Bi}\left(\frac{1}{4\epsilon^{4/3}}\right)}{\left(\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Bi}'\left(\frac{1}{4\epsilon^{4/3}}\right)-\textrm{Ai}'\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Bi}\left(\frac{1}{4\epsilon^{4/3}}\right)\right)},$$ and 
$$c_2 = -\frac{\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)}{\textrm{Bi}\left(\frac{1}{4\epsilon^{4/3}}\right)} $$
Calculating these values of the Airy functions and their derivatives is not necessarily simple.
Substituting back for $W$ gives
$$ W(t) = c_1\textrm{Ai}\left(\frac{1-4\epsilon t}{4\epsilon^{4/3}}\right)+c_2\textrm{Bi}\left(\frac{1-4\epsilon t}{4\epsilon^{4/3}}\right) $$
and so
$$x(t) = \frac{\epsilon^{1/3}e^{-t/(2\epsilon)}\left[-\textrm{Bi}\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Ai}\left(\frac{1-4\epsilon t}{4\epsilon^{4/3}}\right)+\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Bi}\left(\frac{1-4\epsilon t}{4\epsilon^{4/3}}\right)\right]}{\textrm{Ai}\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Bi}'\left(\frac{1}{4\epsilon^{4/3}}\right)-\textrm{Ai}'\left(\frac{1}{4\epsilon^{4/3}}\right)\textrm{Bi}\left(\frac{1}{4\epsilon^{4/3}}\right)}. $$
I do think WKB theory is necessary for the asymptotics (especially if you replace $t$ by $\cos(t)$), and for the equation in terms of $W$ there will be three regions, one where $1/4-\epsilon t>0$, one where it is negative, and a connection region where it is small. The book "Introduction to perturbation methods" by Mark Holmes has a good section on WKB problems with turning points.
