Suppose you have a regular Sturm-Liouville operator
$$
Lf = \frac{1}{w(x)}\left[-\frac{d}{dx}\left(p(x)\frac{df}{dx}\right)+q(x)f\right]
$$
For example, suppose $p \in C^1[a,b]$, $w,q\in C[a,b]$ and $p > 0$, $w > 0$ on $[a,b]$. The domain $\mathcal{D}(L)$ should include some real endpoint conditions such as
$$
\cos\alpha f(a)+\sin\alpha f'(a)=0,\\
\cos\beta f(b)+\sin\beta f'(b) = 0.
$$
(Assume that $\alpha,\beta$ are real so that the endpoint conditions are real.)
By standard existence theorems, there are unique $C^2[a,b]$ solutions of $Lf=\lambda f$ such that
$$
\varphi_{\lambda}(a)=-\sin\alpha,\;\; \varphi_{\lambda}'(a)=\cos\alpha, \\
\psi_{\lambda}(b)=-\sin\beta,\;\;\psi_{\lambda}'(b)=\cos\beta.
$$
The solution $\varphi_{\lambda}$ satisfies the left endpoint condition, while $\psi_{\lambda}$ satisfies the right endpoint condition. These solutions have power series expansions in $\lambda$ because of the fixed endpoint conditions. And you show that the following does not depend on $x$.
$$
w(\lambda)=W_p(\varphi_{\lambda},\psi_{\lambda})(x)=p(x)\{\varphi_{\lambda}(x)\psi_{\lambda}'(x)-\varphi_{\lambda}'(x)\psi_{\lambda}(x)\}
$$
The Wronskian $w(\lambda)$ is an entire function of $\lambda$, and you can show that this function of $\lambda$ vanishes for some $\lambda$ iff $\{ \varphi_{\lambda},\psi_{\lambda} \}$ is a linearly-dependent set of functions, which is equivalent to both functions satisfying both endpoint conditions. That, in turn, is equivalent to the existence of a non-trivial $C^2[a,b]$ solution of $Lf=\lambda f$ which satisfies both endpoint conditions. The zeros of the Wronskian $w$ are the eigenvalues, which cannot cluster because of the identity theorem for holomorphic functions.
If you want to avoid Complex Analysis, you can use Functional Analysis to work in the Banach space $C[a,b]$ and show that $(L-\lambda I)$ has a bounded inverse for $\lambda\notin\mathbb{R}$ that is given by
$$
(L-\lambda I)^{-1}f= \frac{\psi_{\lambda}(x)}{w(\lambda)}\int_{a}^{x}f(x)\varphi_{\lambda}(x)dx
+ \frac{\varphi_{\lambda}(x)}{w(\lambda)}\int_{x}^{b}f(x)\psi_{\lambda}(x)dx.
$$
This is a compact operator because of its smoothing properties, which is shown with a classical application of Arzela-Ascoli; this is done by showing that the integral operator maps bounded sets of functions to equicontinuous sets in $C[a,b]$. So you know a lot about its eigenfunctions and eigenvalues. And you can easily relate the eigenfunctions/eigenvalues of $(L-\lambda I)$ to those of $(L-\lambda I)^{-1}$. You can show that the eigenvalues of $L$ cannot cluster in the finite plane because of properties of the spectrum of compact operators.
You can approach the subject by looking at the zeros of $Lf=\lambda f$ for real $\lambda$. This is another approach, and one of the earliest taken by Sturm.
The differential equation you seek is
$$y''+\lambda y = 0$$
which had the general solution
$$y(x) = A \cos{\sqrt{\lambda} x} + B \sin{\sqrt{\lambda} x}$$
The boundary value given by $y'(0)=y(0)$ implies that $A = B \sqrt{\lambda}$. The boundary value given by $y'(1) + y(1) = 0$ implies that
$$-A \sqrt{\lambda} \sin{\sqrt{\lambda}} + B \sqrt{\lambda} \cos{\sqrt{\lambda}} + A \cos{\sqrt{\lambda}} + B \sin{\sqrt{\lambda}}=0$$
Using the relation derived from the first boundary condition, we get
$$-\lambda \sin{\sqrt{\lambda}} + \sqrt{\lambda} \cos{\sqrt{\lambda}} + \sqrt{\lambda} \cos{\sqrt{\lambda}} + \sin{\sqrt{\lambda}}=0$$
or, after a little algebra,
$$\tan{\sqrt{\lambda}} = \frac{2 \sqrt{\lambda}}{\lambda-1}$$
Here is a plot of the LHS and RHS so you can see the intersections and hence, the square roots of the eigenvalues. Clearly, for the first few, they must be computed numerically. For the $k$th eigenvalue, where $k$ is large,
$$\lambda_k \sim k^2 \pi^2$$
The eigenfunction (unnormalized) corresponding to that eigenvalue is
$$y_k(x) = \sqrt{\lambda_k} \cos{\sqrt{\lambda_k} x} + \sin{\sqrt{\lambda_k} x}$$
The Green's function may then be found from the eigenvalues/eigenfunctions as follows:
$$G(x,x') = \sum_{k=1}^{\infty} \frac{y_k(x) y_k(x')}{\lambda-\lambda_k}$$
where $G$ satisfies
$$\frac{d^2}{dx^2} G(x,x') + \lambda G(x,x') = \delta(x-x')$$
Best Answer
Let $\varphi_{\lambda}(x)=\cos(\sqrt{\lambda}x)$ and $\psi_{\lambda}(x)=\sin(\sqrt{\lambda}(x-1))/\sqrt{\lambda}$. Then $$ \varphi_{\lambda}''=\lambda\varphi_{\lambda},\; \varphi_{\lambda}(0)=1,\; \varphi_{\lambda}'(0)=0 \\ \psi_{\lambda}''=\lambda\psi_{\lambda},\;\psi_{\lambda}(1)=0,\;\psi_{\lambda}'(1)=1. $$ The Green function solution of $(L-\lambda I)g=f$ subject to $g'(0)=0$, $g(1)=0$ is given by $$ g= \frac{\psi_{\lambda}(x)}{w(\lambda)}\int_0^x f(y)\varphi_{\lambda}(y)dy+\frac{\varphi_{\lambda}(x)}{w(\lambda)}\int_x^1f(y)\psi_{\lambda}(y)dy, $$ where $\omega(\lambda)$ is the Wronskian of the solutions $\varphi_{\lambda}$ and $\psi_{\lambda}$ given by $$ \psi_{\lambda}'\varphi_{\lambda}-\varphi_{\lambda}'\varphi_{\lambda} \\=\cos(\sqrt{\lambda}(x-1))\cos(\sqrt{\lambda}x)+\sin(\sqrt{\lambda}x)\sin(\sqrt{\lambda}(x-1)) \\ =\cos(\sqrt{\lambda}(x-1)-\sqrt{\lambda}x)=\cos(\sqrt{\lambda}). $$ There are poles of order $1$ in the resolvent expression at $\sqrt{\lambda_n}=(n+1/2)\pi$ or $\lambda_n=(n+1/2)^2\pi^2$. These are the eigenvalues, and the eigenfunctions are $\psi_{\lambda_n}$.