I want the Green's function of the operator $\frac{d}{dx}\left((1-x^2)\frac{d}{dx}\right)f$, with the boundary conditions $f(0)=f(1)=0$.
I therefore need to solve the equation $\frac{d}{dx}\left((1-x^2)\frac{dG}{dx}\right)=\delta(x-y).$
If $x\neq y$, we have $(1-x^2)\frac{dG}{dx}=c$ with $c$ a constant and $$ G=\frac{c}{2}\int\frac{dx}{1-x}+\frac{c}{2}\int\frac{dx}{1+x}=\frac{c}{2}\log\left(\frac{1+x}{1-x}\right)+d.$$
This form is valid both for $x<y$ and $x>y$. Wecan use $c_1,d_1$ in the first domain and $c_2,d_2$ in the other.
Now, I muts impose the boundary conditions. Vanishing at $x=0$ implies $d_1=0$. But vanishing at $x=1$ implies $c_2=d_2=0$. So $G(x,y)=\frac{c_1}{2}\log\left(\frac{1+x}{1-x}\right)$ if $x<y$ and $G(x,y)=0$ if $x>y$.
If I impose continuity at $x=y$, I get $c_1=0$. So $G=0$ everywhere.
This can't be right. What gives?
You cannot impose traditional endpoint conditions for this equation because the equation is singular at $x=\pm 1$. The solution of $((1-x^2)f')'=0$ have the form that you discovered, which is $$ f = C\frac{1}{2}\ln\frac{1+x}{1-x}+D. $$ So there are two arbitrary constants associated with $f$ at each of the two endpoints, but imposing $f(0)=0$ (or $f(1)=0$) requires $C=D=0$. This is because the ODE is singular at $x=\pm 1$. The above may be treated as an asymptotic for homogeneous solutions near the endpoints, and the constants take the place of $f(0),f'(0),f(1),f'(1)$ for a non-singular equation. Physicists typically require boundedness near the endpoints, which forces $C=0$ and allows $D$ to vary, which allows you to construct the Green function.
(As a picky note, it is "Green function" and not "Green's function". It sounds wrong because green is a color, but we say Bessel function or Legendre function, etc., and not Bessel's function(s), etc.)