I'm attempting to numerically solve this differential equation using MATLAB's ode45, presuming a stable solution, i.e., $y$ remains finite as $x\to+\infty$. ($g(x)$ is a pending function determined by my experimental data. The value of $g(x)$ is always real and positive.)
Here is my confusion:
If $g(x)$ is just a constant $g_0$, the analytical solution is, of course, $y=C_1\exp(\sqrt{g_0}x)+C_2\exp(-\sqrt{g_0}x)$. ($C_1, C_2$ are constants.) We can also obtain a stable numerical solution by setting the initial values as $y(0)=1; y'(0)=-\sqrt{g_0}$.
If $g(x)$ varies slowly with $x$, analytically the equation can still be solved using the WKB method. The zeroth-order solution is $y\simeq C_1\exp\left(\int_0^x\sqrt{g(x')}dx'\right)+C_2\exp\left(-\int_0^x\sqrt{g(x')}dx'\right)$. However, for the numerical approach, I could no longer get a stable solution by merely changing the initial values in the ode45 solver. I've tried $y(0)=1; y'(0)=-\sqrt{g(0)}$ but obtained an unstable blue line in the figure below. The result also does not match the WKB solution $\exp\left(-\int_0^x\sqrt{g(x')}dx'\right)$ indicated by the red dashed line.
Is there any way to avoid the instability driven by the $+\sqrt{g(x)}$ exponential branch? This question has puzzled me so much, and any valuable advice will be appreciated. (I'm also wondering whether this question is more appropriate to be posted here or on the Computational Science Stack Exchange.)