Function $y(x)$ and $g(x)$, defined on interval $[0,1]$, are related by linear differential equation that can be read as a second order or as a first order ODE
\begin{equation} g~y’’+\frac{1}{2}~g’~y’+\frac{1}{4}~((g-1)/x)’~y~=~0, \tag{1} \end{equation} with boundary conditions \begin{equation} y(1)=y_1,~~~y’(1)=y’_1,~~~g(1)=g_1. \tag{2} \end{equation} I have found in the literature following method of obtaining solution - one postulates some function, say $\Phi(x)$, and calls the function \begin{equation} y=a~\Phi(b~x), \tag{3} \end{equation} with \begin{equation} a\equiv\frac{y_1}{\Phi(b)},~~~b\equiv \frac{\Phi(b)}{\Phi'(b)}\frac{y'_1}{y_1}. \tag{4} \end{equation} solution of equation $(1)$ to corresponding boundary conditions.
Using it one solves again equation $(1)$ this time for the function $g$ and states that the postulated function $y$ is the solution of the second order equation $(1)$ with variable coefficients defined by the function $g$.
I argue that it is not correct because the solution space of the second order differential equation is a two dimensional vector space with base built by two linearly independent particular solutions . The true solution has to be combination of them, say functions $\Phi(x)$ and $\Psi(x)$, where the second is derived from $\Phi(x)$ by the method of reduction order, both satisfying equation $(1)$, and the coefficients of the linear combination have to be determined by the boundary conditions $(2)$. For me the solution looks like \begin{equation} y=c_1~\Phi(x)+c_2~\Psi(x), \tag{5} \end{equation}
with \begin{equation} c_1~=~\frac{y_1~\Psi'_1-y'_1~\Psi_1}{\Phi_1~\Psi'_1-\Phi'_1~\Psi_1},~~~c_2~=~-\frac{y_1~\Phi'_1-y'_1~\Phi_1}{\Phi_1~\Psi'_1-\Phi'_1~\Psi_1}. \tag{6} \end{equation} I am not mathematician and the physicists that have used the method I claim to be wrong are numerous and well-known - therefore my doubts about it.