Sign of solution to (in)homogeneous linear ODE Let $N \geq 3$ be a positive integer and $A >0, B \geq 0$ be two constants. Let $y: (0,\infty) \to \mathbf{R}$ be a solution to the following linear, inhomogeneous ODE: $y''(x) + \frac{N-1}{x} y'(x) + (\frac{N-1}{x^2} + A) y(x) = - B$ for all $x \in (0,\infty)$ with initial values $y(0) = 0, y'(0) = 0$. (In the case $B = 0$ the equation is homogeneous).
Question: What is the sign of $y$ near $x = 0$? Does it have a sign regardless of $A,B$? What if $A = B/2$, or $B = 0$?
I would also be grateful for comments on a special case or a pointer to a chapter in the literature where this is treated.
Attempts. The ODE has a regular singular point at $x = 0$, which makes it amenable to a solution via formal power series. After consulting computer programs, it seems that the solutions of the homogeneous ODE are $x^{-\tau} J_{\sqrt{k}/2}$ and $x^{-\tau} Y_{\sqrt{k}/2}$ where $\tau > 0$, $k \in \mathbf{Z}_{>0}$ and $J_{\sqrt{k}/2},Y_{\sqrt{k}/2}$ are Bessel functions. I had trouble recovering this expression; moreover when $B \neq 0$ it seems too complicated to recover the sign.
On the other hand, after multiplying the equation by $x^{N-1}$ one obtains an equation in Sturm--Liouville form $(x^{N-1} y')' + x^{N-1}(\frac{N-1}{x^2} + A)y = - Bx^{N-1}$, or after flipping the sign $-(py')' + qy = Bx^{N-1}$ with $p = x^{N-1}$ and $q = -x^{N-1}(\frac{N-1}{x^2} + A)$. As far as I understand Sturm--Liouville theory ought to give some information about the sign of the solutions, but I was not able to find an answer, even in the homogeneous case. (Presumably if the solution to the homogeneous ODE were negative, then by comparison the same ought to be true when $B > 0$?)
 A: Your equation has a particular entire solution of the form
$$y(x)=\sum_{n=1}^\infty c_nx^{2n},\quad c_1=-B/(3N-1).$$
This solution can be obtained by substituting this series to the equation
and determining all coefficients one-by-one. A solution of this form is unique,
and it satisfies your initial conditions. All coefficients $c_k$ can be explicitly determined from a simple recurrent relation.
It is negative near $0$.
The general solution is obtained as a sum of this one and a linear combination
of solutions of the homogeneous equation. Adding a linear combination of the homogeneous equation may not change your initial condition if the real part of the exponent at $0$
of this combination is $>1$. If such a solution exists, then the solution of your initial value problem is not unique. Let us determine when it exists. The characteristic equation is
$$\rho^2+(N-2)\rho+N-1=0.$$
Its larger solution is $(2-N+\sqrt{N^2-8N-8)})/2$. This cannot be $\geq 1$, and if it
is non-real, then real part is negative. So a solution of the homogeneous equation
cannot affect the sign of $y(x)$ near $0$.
The conclusion is that your initial value problem always has a solution,
(for your range of parameters) and this solution is negative near $0$;
even in the case when the solution of your initial values problem is not unique, all solutions have the same sign near $0$.
