I am not at all an expert in Sturm-Liouville theory, but I ended up on the following Singular Sturm Liouville problem: \begin{equation}\label{1} (1) \ \ \ \ \ \ \ \ \ \ \ y''(t)+\frac{\theta'(t)}{\theta(t)}y'(t)+\lambda y(t)=0\,, \ \ \ t\in(0,R) \end{equation} which can be rewritten as $$ (\theta(t) y'(t) )'+\lambda\theta(t)y(t) =0\,, \ \ \ t\in(0,R). $$ The weight function is strictly positive and smooth on $(0,R)$, and vanishes at $t=0$ and $t=R$, which are then singular endpoints. I don't have much information on $\theta$ except that it has polynomial behavior near the endpoints, namely \begin{equation}\label{2} (2) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \theta(t)\sim_{t\rightarrow 0^+} t^a\ \ \ {\rm and}\ \ \ \theta(t)\sim_{t\rightarrow R^-} (R-t)^a. \end{equation} The crucial point is that I need the boundary conditions (which I can't impose in principle) $$ y'(0)=y'(R)=0 $$ In specific cases I know how to handle this (e.g., when $\theta(r)=\sin(r\pi/R)^a$ by suitable change of variables I end up on a Jacobi equation).
In the general case, from Sturm-Liouville theory (if I am not wrong, otherwise please correct me), if $1\leq a<3$, the endpoints are of LC type, while they are LP type if $a\geq 3$.
In the LC case, the natural choice for boundary conditions for me is $\theta(t)y'(t)\rightarrow 0$ as $t\rightarrow 0^+$ and $R^-$.
In the LP case such boundary conditions are in a sense authomatic and no other options are available. In both cases I have discrete spectrum, made of simple non-negative eigenvalues $\lambda_n$ diverging to $+\infty$ (due to the behavior of $\theta$ near the endpoints I can see that in the LP case the essential spectrum is empty).
The only point that is left to prove is: do the eigenfunctions $y_n$ satisfy also $y_n'(0)=y_n'(R)=0$ (and not only $\theta(t)y_n'(r)\rightarrow 0$ as $t\rightarrow 0^+$ and $R^-$)?
An example to motivate this (perhaps experts have already the answer): in the simple case $$ y''(t)+\frac{a}{t}y'(t)+\lambda y(t)=0\,,\ \ \ t\in(0,1) $$ where $t=0$ is the only singular endpoint, I can put the boundary condition $t^ay'(t)\rightarrow 0$ as $t\rightarrow 0$ (and any "classical" BC at t=1) if $1\leq a<3$. The same b.c. at $t=0$ is authomatic if $a\geq 3$. Then I have in both cases that $y_n(t)=x^{\frac{1-a}{2}}J_{\frac{a-1}{2}}(\sqrt{\lambda_n}x)$, where $\lambda_n$ are provided by imposing the classical b.c. at $t=1$ (for $\lambda_0=0$, $y_0=1$).
These eigenfunctions satisfy the condition $t^ay_n'(t)\rightarrow 0$ (which I have imposed in the case $1\leq a<3$, and which is authomatic for $a\geq 3$). Not only, they satisfy $y_n'(0)=0$ (but this I can recover from well-known theory of Bessel functions, not from the form of the problem).
To summarize, my question is the following:
"Do the eigenfunctions $y_n$ of (1) with weights satisfying (2), and with $\theta(t)y_n'(t)\rightarrow 0$ as $t\rightarrow 0^+$ and $R^-$ satisfy also the condition $y_n'(0)=y_n'(R)=0$?"
This seems quite reasonable and intuitive to me. I would appreciate some reference where to find this kind of result.