Consider the set of integrals $$I_n=\int_0^{F(\varphi,m)}(c+\operatorname{sn}u)^{-n}\,du=\int_0^{\sin\varphi}\frac1{(c+t)^n\sqrt{(1-t^2)(1-mt^2)}}\,dt$$ This is the series denoted $\sigma_n$ in Byrd and Friedman's Handbook of Elliptic Integrals for Engineers and Physicists. They give formulas that suffice to compute $I_n$ for $n\le2$ (330.00–04, .50, .51).
However, in my work I have come across the need to evaluate $I_3,I_4$ and higher indices of $I_n$. I cannot see a way to convert these integrals into other integrals covered in B&F's "reduction of integrals of elliptic functions" section, and I have a feeling that I will need another recurrence to handle them. How can I derive a recurrence for $I_n$ that allows me to solve $I_3,I_4$ and so on, or can I solve it by any other means? More generally, how did Byrd and Friedman come up with all those recurrences in the first place?
Using the standard reduction formula (B&F 250.02) applied to the algebraic form of $I_n$ gives for all $n\ge2$ $$I_n=\frac1{(n-1)(1-c^2)(1-mc^2)}\left(\frac1{c^{n-1}}-\frac{\operatorname{cn}u\operatorname{dn}u}{(c+\operatorname{sn}u)^{n-1}}+\sum_{j=1}^4C_jI_{n-j}\right)$$ where $$\begin{align} C_1&=c(2n-3)(2c^2m-m-1)\\ C_2&=-(n-2)(6c^2m-m-1)\\ C_3&=2cm(2n-5)\\ C_4&=-m(n-3) \end{align}$$ provided that $c^2\ne1$ and $c^2\ne\frac1m$. Otherwise, B&F 250.03 applies: $$I_n=\frac1{c(2n-1)(2c^2m-m-1)}\left(\frac{\operatorname{cn}u\operatorname{dn}u}{(c+\operatorname{sn}u)^n}-\frac1{c^n}+\sum_{j=1}^3D_jI_{n-j}\right)$$ where $$\begin{align} D_1&=(n-1)(6c^2m-m-1)\\ D_2&=-2cm(2n-3)\\ D_3&=m(n-2) \end{align}$$ It is perhaps because of the need to distinguish between $P(c)\ne0$ and $P(c)=0$ (where $P$ is the polynomial under root), as well as the complicated factors involved, that B&F does not list recurrence relations for the $I_n$ here.