In the theory of ellipsoidal harmonics, Lame functions of the second kind $F_n$ arise as the second linearly independent solution (the first being Lame functions of the first kind $E_n$) of Laplace's equation in ellipsoidal coordinates which is the so called Lame's equation. Lame functions of the second kind contain the improper elliptic integral $I_n(\rho)$ and for this reason they are obtainable through numerical approximations e.g. Gauss-Legrendre quadrature. A really comprehensive book on the theory of ellipsoidal harmonics is https://books.google.gr/books?isbn=0521113091  .

For "investigational" purposes my goal is to find an expansion in powers of 1/ρ (and its first 2 or 3 terms) of the elliptic integral:
\begin{equation}I_n(\rho)=\int_\rho^{+\infty}\frac{dt}{(E_n(t))^2\sqrt{t^2-h_2^2}\sqrt{t^2-h_3^2}},\quad \rho \ge h_2\end{equation}
where $E_n(t)$ is a polynomial  in descending powers of n  of the form:
\begin{equation}
E_n(t)=\sum_{k=0}^\infty a_kt^{n-2k}, \quad a_0\neq 0\end{equation}
where for negative powers summation is terminated.
A first step to solution is to set $t\rightarrow \frac{h_2}{t}$ and this leads to
\begin{equation}
    I_n(\rho)=\int_0^\frac{h_2}{\rho}\frac{t^{2n}/h_2^{2n}}{(E_n(t))^2\sqrt{1-t^2}\sqrt{1-\frac{h_3^2}{h_2^2}t^2}}dt. 
\end{equation}
Using binomial theorem
\begin{align}
(1-t^2)^{-1/2}&=\sum_{n=0}^{+\infty}(-1)^n\binom{-\frac{1}{2}}{n}t^{2n},\\
(1-\frac{h_3^2}{h_2^2}t^2)^{-1/2}&=\sum_{n=0}^{+\infty}(-1)^n\binom{-\frac{1}{2}}{n}\frac{h_3^{2n}}{h_2^{2n}}t^{2n},
\end{align}
and multiplying gives the series expansion 
\begin{equation}\sum_{n=0}^{+\infty}\sum_{k=0}^n(-1)^n\binom{-\frac{1}{2}}{n-k}\binom{-\frac{1}{2}}{k}\frac{h_3^{2k}}{h_2^{2k+1}}t^{2n}\end{equation} so:
\begin{equation}I_n(\rho)=\int_0^{h_2/\rho}\left[\frac{t^{2n}/h_2^{2n}}{(E_n(t))^2 }\sum_{n=0}^{+\infty}\sum_{k=0}^n(-1)^n\binom{-\frac{1}{2}}{n-k}\binom{-\frac{1}{2}}{k}\frac{h_3^{2k}}{h_2^{2k+1}}t^{2n}\right] dt .\end{equation}

And here starts the problem:

The reciprocal of the square of $E_n$ is:
\begin{equation}
\frac{1}{E_n[(t)]^2}=\sum_{n=0}^{+\infty}c_nt^n
\end{equation}
where:
\begin{align}
c_0&=\frac{1}{a_0}\\
c_n&=-\frac{1}{na_0^2}\sum_{k=0}^{n-1}(2n-k)a_{n-k}c_k, \quad n\ge 1
\end{align}

If everything up to now is correct then how do i proceed? 
Hope you like the subject and offer some help on the matter. 
Thank you very much :)