I have a system of two linear differential equations with non-constant coefficients (which involve hypergeometric functions),
$$ u_1'(y) = \alpha_1(y) u_1(y) + \beta_1(y) u_2(y) \\ u_2'(y) = \alpha_2(y) u_1(y) + \beta_2(y) u_2(y) $$
that I want to integrate numerically in the interval $y\in[0,1]$. The coefficients are singular at the boundary of the interval (i.e., at $y=0$ and $y=1$).
As we have a linear system of 2 equations, the solutions conform a 2-dimensional vector space of functions, and a general solution can be written as
$$ u(y) = C_\text{I} u^\text{I}(y)+C_\text{II}u^\text{II}(y) $$
where $u=(u_1,u_2)$. By construction, every particular solution (that is, every choice of $C_\text{I}$ and $C_\text{II}$) vanishes when $y\to1$, but there's only a 1-dimensional subspace of solutions that gives finite value for $u_1$ and $u_2$ when $y\to0$.
I have the asymptotic expansion of the solution around $y=0$ and $y=1$. These series give us $u_1$ and $u_2$ around $y=0$ depending on some constants $C_\text{I}^{(0)}$ and $C_\text{II}^{(0)}$, and around $y=1$ with another set of constants $C_\text{I}^{(1)}$ and $C_\text{II}^{(1)}$.
I know a relation between $C_\text{I}^{(0)}$ and $C_\text{II}^{(0)}$, $G(C_\text{I}^{(0)}$,$C_\text{II}^{(0)})=0$, such that the solution doesn't explode on $y=0$, but I want to know the equivalent choice of $C_\text{I}^{(1)}$ and $C_\text{II}^{(1)}$ for the expansion around $y=1$, that allows me to integrate the system from $y=1$, ensuring that the solution will be finte at $y=0$.
Any insight on how could I achieve this in Mathematica? I have been looking about the Shooting Method, so I could shot solutions from $y=1$ and pick the one (up to a multiplicative constant factor) that is finite at $y=0$, but I am not sure on how could I implement this.
