In order to clarify my problem, I start with an example $$ y'-\frac{1}{2}\left(\frac{1}{x}+\frac{1}{x-1}\right)y=0 $$ which has two singular points, $0$ and $1$. The exact solution is not difficult to be obtained, but I would like to do the local analysis around each singular points, and find the relations between them.
Since the coefficient of the differential equation is analytical in $0<|x|<1$ and $0<|x-1|<1$, the solutions around each singular points can be represented by $$ y_0= x^{\alpha} f(x), \qquad y_1=(x-1)^{\beta} g(x) $$ where $f$ and $g$ are single-valued functions in corresponding ranges.
To perform the local analysis, I am going to use the seise expansions for $f$ and $g$, i.e. $$ f=\sum_{n=0} a_n x^n,\qquad g=\sum_{n=0} b_n x^n $$ Then by substituting the ansatz into the original equation, I can fix the $\alpha$ and $\beta$, the coefficients of the series and the relation between the solutions at each points $$ a_n = (-1)^n b_n $$
My question: For a given differential equation $$ y'+p(x) y=0 $$ where $p(x)$ has to two simple poles, say $x_0$ and $x_1$, if I know an asymptotic (series) solution around $x_0$, can we obtain the solution around $x_1$ without solving the equation?
Step one: if $p(x)$ is an entire function, thus around some regular point $x_0$, one has a series solution of the differential equation $$ y(x)={\bf A}\cdot {\bf X}, $$ where ${\bf A}=\{a_0, a_1, \ldots\}$ and ${\bf X}=\{1, (x-x_0), \ldots\}$ are two infinite dimensional vectors. The radius of convergence is $R$. For another regular point $x_1$ located in the convergence range of above series, we can obtain a series solution at $x_1$ directly by $$ y(x)={\bf A}\cdot {\bf J} \cdot{\bf Y}, $$ where ${\bf Y}=\{1, (x-x_1), \ldots\}$ and the component of the matrix ${\bf J}$ are $$ J_{ij}=\binom{i}{j}\Delta^{i-j}. $$ Here $\Delta=|x_1-x_2|<R$. If $x_1$ is not located in the convergence range, one can make a chain, such that $$ y(x)={\bf A}\cdot J_1 \cdot\ldots \cdot J_n \cdot {\bf Y}. $$
Step two: if $p(x)$ has two simple poles, the space is extended, at $x_0$ one has $$ y(x) = \left( \tilde A \; A \right) \cdot \begin{pmatrix} \tilde X \\ X \\ \end{pmatrix} $$ where $\tilde A=\{0,\ldots,0,a_{-1}\}$ and $\tilde X=\{(x-x_{0})^{-\infty},\ldots,(x-x_{0})^{-1}\}$, then around $x_1$, we have $$ y(x) = \left( \tilde A \; A \right) \cdot \begin{pmatrix} {\bf K} & {\bf 0}\\ {\bf 0} & {\bf J} \end{pmatrix} \cdot \begin{pmatrix} \tilde Y \\ Y \\ \end{pmatrix} $$ where the component of ${\bf K}$ is $$ K_{ij}= \delta_{-1,i} \Delta^{j-1}. $$ Here $\delta_{-1,i}$ is Kronecker delta. As a result, one has $$ y(x) = \left( \tilde A \cdot {\bf K}\;\; A \cdot {\bf J} \right) \cdot \begin{pmatrix} \tilde Y \\ Y \\ \end{pmatrix} $$