I am numerically solving the following boundary value problem:
$$ \begin{align} y'_1 &= F(x)+x^2\rho(x) y_2\\ y'_2 &= -c\frac{\rho(x)}{x^2 P(x)}y_1 \end{align} $$ with the boundary conditions $y_1(x=0)=y_1(x=1)=0$.
where $x$ is the independent variable, $\rho, P, F$ are functions of $x$, and $c$ is a positive constant.
I have solved this system of equations for a large number of combinations of $\rho, P, F$. I have come across a curious situation in which a tiny change to the independent variable $x$ results in a change of sign to the solutions $y_1,y_2$.
I am wondering if anyone has advice as to how I can go about understanding this apparent instability (and potentially remove it)?