I have
$$ \partial_{\epsilon}^{2}g+\frac{R^{4}\left(R^{2}+3\epsilon^{2}\right)}{4\epsilon^{3}\left(R^{2}-\epsilon^{2}\right)}\partial_{\epsilon}g+\frac{R^{4}\left(5R^{2}+3\epsilon^{2}\right)}{4\epsilon^{2}\left(R^{2}-\epsilon^{2}\right)^{2}}g=0 $$
and I've got one of the solutions,
$$ g_{1}=\frac{1}{\pi R^{2}}\left(2\left(R^{2}-\epsilon^{2}\right)E\left(\frac{\epsilon}{R}\right)\right) $$
where $E$ is the complete elliptic integral of the second kind.
I've tried to use the method of reduction of order (http://tutorial.math.lamar.edu/Classes/DE/ReductionofOrder.aspx) to find the other solution, though it seems that the $v(x)$ term (where the new solution has the form $v(x)g_{1}$) in the new DE doesn't drop out (as it should in this method).
Can anyone help me find the other solution --OR-- show that there is only one unique solution ($= g_{1}$)?
Note: If you want to know how I got the first solution, look at the relation of the elliptic integral to the hypergeometric differential equation, working backwards....
we have
$$ g^{\star}(\epsilon)=\frac{1}{\pi R^{2}}\left(2\left(R^{2}-\epsilon^{2}\right)E\left(\frac{\epsilon}{R}\right)\right) $$
We notice that the complete elliptic integral of the second kind, E, is related to Gauss' hypergeometric function $_{2}F_{1}$ (which we write as F) $$F(a,b;c;z)=\sum_{n=0}^{\infty}\frac{\left(a\right)_{n}\left(b\right)_{n}z^{n}}{\left(c\right)_{n}n!}$$
where$$\left(\alpha\right)_{\gamma}\equiv\frac{\Gamma\left(\alpha+\gamma\right)}{\Gamma\left(\alpha\right)}$$
is the Pochhammer symbol. We also have
$$E\left(k\right)=\frac{\pi}{2}F\left(-\frac{1}{2},\frac{1}{2};1;k^{2}\right)$$
such that
$$F\left(-\frac{1}{2},\frac{1}{2};1;\frac{\epsilon^{2}}{R^{2}}\right)=\frac{R^{2}g^{\star}(\epsilon)}{R^{2}-\epsilon^{2}}$$
and since we have the hypergeometric differential equation
$$z(1-z)\partial_{z}^{2}w+\left(c-\left(a+b+1\right)z\right)\partial_{z}w-abw=0$$
whenever
$w=F\left(a,b;c;z\right)$, we have
$$ \frac{\epsilon^{2}}{R^{2}}(1-\frac{\epsilon^{2}}{R^{2}})\partial_{\frac{\epsilon^{2}}{R^{2}}}^{2}F+\left(1-\frac{\epsilon^{2}}{R^{2}}\right)\partial_{\frac{\epsilon^{2}}{R^{2}}}F+\frac{F}{4}=0 $$
which simplifies to a nonlinear second order differential equation in $g^{\star}$
$$\partial_{\epsilon}^{2}g^{\star}+\frac{R^{4}\left(R^{2}+3\epsilon^{2}\right)}{4\epsilon^{3}\left(R^{2}-\epsilon^{2}\right)}\partial_{\epsilon}g^{\star}+\frac{R^{4}\left(5R^{2}+3\epsilon^{2}\right)}{4\epsilon^{2}\left(R^{2}-\epsilon^{2}\right)^{2}}g^{\star}=0$$
First of all, this is an ODE, not a PDE, because you're differentiating with respect to just one independent variable.
Second, I don't get $g_1$ as a solution. This is the most likely reason why the $v(x)$ term doesn't drop out of the reduction of order. Reduction of order should always find a second solution if your first solution is correct and nonzero.
EDIT: In fact, what I find is that $g_1$ is a solution of a slightly different equation:
$$ g''(\epsilon) + \dfrac{R^2 + 3 \epsilon^2}{\epsilon (R^2 - \epsilon^2)} g'(\epsilon) + \dfrac{5 R^2 + 3 \epsilon^2}{(R^2 - \epsilon^2)^2} g = 0 $$ This is using Maple's parametrization of the elliptic integral; Mathematica's is somewhat different.
Third, since this is a second order linear DE, there are always two linearly independent solutions.
Maple 18, by the way, gets solutions in terms of the HeunC function.