This may be a naive question and I apologize if it's not stated concisely and rigorously. I try to pose it first, and then explain my thoughts and comments on it; so please bear with me.
Let $\lambda_n$ be the eigenvalues of the Sturm-Liouville problem $Ly_n=\lambda_n y_n$ with homogeneous boundary conditions $y_n(x)=0$ at $x=0$ and $x=1$. Now consider $Ly_\lambda=\lambda y_\lambda$ with inhomogeneous boundary conditions $y_\lambda(0)=0$ and $y_\lambda(1)=1$. The questions is: Under what conditions is the following statement true?
$y_\lambda$, as a function of $\lambda$, has poles at $\lambda_n$ [for almost all $x\in(0,1)$].
This statement is true for $L=-d^2/dx^2$ where $\lambda_n=n^2\pi^2$ and $y_\lambda = \sin\sqrt\lambda x / \sin\sqrt\lambda$. I can argue for it in a suggestive way that may be generalizable to other cases. The general solution of $Ly=\lambda y$ for $\lambda\neq\lambda_n$ is $y(0) y_{10}(x) + y(1) y_{01}(x)$, where $y_{01} = \sin\sqrt\lambda x/ \sin\sqrt\lambda$ satisfies the inhomogeneous boundary condition $y_{01}(0)=0$ and $y_{01}(1)=1$. Similarly $y_{10} = \cos\sqrt\lambda x - \cos\sqrt\lambda y_{01}$ satisfies $y_{10}(0)=1$ and $y_{10}(1)=0$. Now I expect to recover $y_n(x)$ from either $y(0) y_{10}(x)$ or $y(1) y_{01}(x)$ as $\lambda\to\lambda_n$, if I properly carry out a limiting procedure to approach the homogeneous boundary conditions by sending $y(0)$ or $y(1)$ to zero. For example, in the current case, I can take $y(1) \equiv c(\lambda) = \sin\sqrt\lambda$. More generally, and adding an explicit label to denote the $\lambda$-dependence, I expect $$\lim_{\lambda\to\lambda_n} c(\lambda) y^{(\lambda)}_{01}(x) = y_n(x).$$ [This $y^{(\lambda)}_{01}$ is the same as the $y_\lambda$ I defined above.] Since $y_n$ is almost always nonzero and I have $\lim_{\lambda\to\lambda_n} c(\lambda) = 0$, I get $\lim_{\lambda\to\lambda_n} y^{(\lambda)}_{01} = \infty$, as claimed. Of course, this was by no means a rigorous argument.
I have a couple of additional comments. First, being a physicist, I grant all nice smoothness assumptions on various functions. Second, I am also interested in situations where one of the boundary conditions is of Neumann type (on the derivative of $y$). Last, but not least, I'd like to know if and when a similar statement holds for a PDE on a domain $\Omega$. That is, for a general (maybe elliptic?) linear operator $L$, let $y_n$ be the solutions of $Ly_n=\lambda_n y_n$ with a homogeneous boundary condition $y(\partial\Omega)=0$ (or a Neumann boundary condition), where $\partial\Omega$ has two disconnected components $\partial\Omega_0$ and $\partial\Omega_1$. Then let $y_\lambda$ be the solution of $Ly_\lambda=\lambda y_\lambda$ with $y_\lambda(\partial\Omega_0)=0$ and $y_\lambda(\partial\Omega_1)=1$. How are the singularities of $y_\lambda$ related to the spectrum of $L$?
If $u$ is any non-trivial solution of $Lu=\lambda u$ for which $u(0,\lambda)=0$, then $$ v(x,\lambda)= \frac{u(x,\lambda)}{u(1,\lambda)} $$ is a solution of $Lv=\lambda v$ for which $v(0,\lambda)=0$, $v(1,\lambda)=1$, unless $u(1,\lambda)=0$, in which case $v(1,\lambda)$ has a pole of order $1$ at $\lambda=1$.