Asymptotic expansion of inhomogenous differential equation Consider
$$\tag{1}
y'(x)+y(x)=\frac{1}{x}
$$
For reference, the exact solution is
$$\tag{2}
y(x)=e^{-x}(C+\operatorname{Ei}(x))
$$
Where $\operatorname{Ei}$ is the exponential integral and $C$ is the integration constant. I want to study the/a particular solution of (1) as $x \to 0^+$. Using dominant balance, I have found $y \sim \ln x$. This matches the logarithmic singularity carried by $\operatorname{Ei}$.
Question: is it possible to say anything about the next to leading order terms of the particular solution using asymptotic analysis? Ie. by manipulating (1) and not just expanding (2). Expanding $\operatorname{Ei}$ suggests the next term should be $\gamma$, Euler's constant.
My thoughts: I think the answer is 'no' because a constant term could be absorbed into the constant of integration, $C$. This leaves me vaguely uneasy.
Working:
There are three dominant balances to consider  in (1). The consistent one uses $y\ll x^{-1}$
$$
y'\sim x^{-1} \qquad, \qquad x \to 0
$$
This can be directly integrated
$$
y(x)=\ln(x)+A(x)
$$
Where $A'(x)\ll x^{-1}$ as $x\to 0$. We may substitute into (1) to find a differential equation for $A$
$$
A'+A=-\ln(x)
$$
The dominant balance that neglects $A$ leads to
$$
A(x)=-x \ln(x)+x+B(x)
$$
Where $B'\ll \ln x$. The differential equation for $B$ is
$$\tag{3}
B'+B=x \ln(x)-x
$$
And now there are two consistent dominant balances to consider, and the RHS vanishes at zero. At this point I think something must have been missed, because no constant terms appear in our expansion (4) between the singular $\ln(x)$ and finite $x \ln(x)$ terms
$$\tag{4}
y \sim \ln(x) - x \ln(x) +x \qquad ,\qquad x \to 0
$$
After playing around with it, I noticed that continuing to neglect the non-derivative terms in (3) and beyond leads to repeated integrals over $x \ln(x)$, which can be done. Spotting the pattern then summing up the terms, I find
$$
y \sim p(x)+e^{-x}\ln(x) \qquad , \qquad x \to 0
$$
Where $p(x)$ is a series in only positive powers of $x$.
Context: The equation (1) comes from this question.
 A: I'm not 100% sure what you are looking for but let me try something. First of all, identifying the constant term in the asymptotic expansion without giving an initial data seems impossible. The reason is because the constant depends on the initial data and is not a universal feature of the ODE itself. I shall briefly touch upon this matter at the end of the answer.

In this answer, let us consider the equation
$$
y'+y=f
$$
for a general function $f(x)$ so as to give a feeling that the approach described below does not require an explicit solution (we'll come back to the case of $f(x)=1/x$ later on); but we require the knowledge of the asymptotic structure of $f(x)$ near $x=0$: we suppose that an asymptotic expansion of the form
$$
f(x)\sim \frac{a_{-1}}{x}+\sum_{n=0}^{\infty}a_n x^n \quad (x\to 0)
$$
holds. The expansion could have started from $a_{-m}x^{-m}$ with $m\geq 2$, but we take $m=1$ for brevity.
Now, looking at the homogeneous terms in the equation $y'+y=f$, it is natural to introduce a new dependent variable $u$ by
$$
y=ue^{-x}.
$$
This is just the method of variation of constants as mentioned in the comment by TheSimpliFire. Then the equation for $u$ is
$$
u'=e^x f.
$$
Let us set $g(x)=e^{x}f(x)$. Since the Taylor expansion of $e^x$ is well-known, we can calculate the coefficients $b_n$ appearing in the following asymptotic expansion of $g(x)$:
$$
g(x)\sim \frac{b_{-1}}{x}+\sum_{n=0}^{\infty}b_n x^n \quad (x\to 0).
$$
For example, $b_{-1}=a_{-1}$ and $b_0=a_{-1}+a_0$.
Next, integrating the equation $u'=g$ from $x$ to $1$, we get
$$
u(x)=u(1)+\int_{1}^{x}g(x')\, dx'.
$$
Since integrations preserve asymptotic expansions (as opposed to derivatives), we can continue the calculation as follows:
\begin{align}
u(x)
& =u(1)+\int_{1}^{x}\frac{b_{-1}}{x'}\, dx'+\int_{0}^{x}\left( g(x')-\frac{b_{-1}}{x'} \right)\, dx'-\int_{0}^{1}\left( g(x')-\frac{b_{-1}}{x'} \right)\, dx' \\
& \sim b_{-1}\log x+\text{const.}+\sum_{n=0}^{\infty}\frac{b_n}{n+1}x^{n+1}.
\end{align}
In this way, we can obtain an asymptotic expansion of $u(x)$, and of course, of $y(x)=e^{-x}u(x)$.
Finally coming back to the case of $f(x)=1/x$, we get
$$
b_{-1}=1, \quad b_n=\frac{1}{(n+1)!} \quad (n\geq 0).
$$
So our conclusion is
$$
y(x)\sim e^{-x}\left( \ln x+\text{const.}+\sum_{n=1}^{\infty}\frac{x^n}{n\cdot n!} \right).
$$

Let us finally comment on the constant term in the asymptotic expansion. It can be calculated using the knowledge of $u(1)$. In fact, we have
$$
\text{const.}=u(1)-\int_{0}^{1}\left( g(x')-\frac{b_{-1}}{x'} \right)\, dx'.
$$
If the asymptotic expansion of $g(x)$ is converging (which we have not assumed so far), we have
$$
\text{const.}=u(1)-\sum_{n=0}^{\infty}\frac{b_n}{n+1}x^{n+1}.
$$
In the case of $f(x)=1/x$, we would get
$$
\text{const.}=u(1)-\operatorname{Ei}(1)+\gamma,
$$
which is consistent with the exact solution in your question.
