I keep seeing it claimed that the general second-order Fuchsian equation with four singularities can be transformed to the Heun equation, but I have never seen anyone explicitly write out the steps, with only a brief mention that Möbius transformations are involved to send three of the singularities at arbitrary position to the standard $0, 1, \infty$.
I am more familiar with Riemann $P$ and how to reduce it so that the solution can be expressed in terms of the Gauss hypergeometric function ${}_2 F_1$, so I thought I'd try to see if what works there can be adapted so that I have a solution in terms of the "local" Heun function, $\operatorname{\mathit{H\ell}}\left(a,q;\alpha,\beta,\gamma,\delta;z\right)$.
As I understand it, the Fuchsian equation with four singularities should look like the following:
$$\frac{\mathrm d^2 w}{\mathrm dz^2}+\left(\frac{1-a_1-a_2}{z-\alpha}+\frac{1-b_1-b_2}{z-\beta}+\frac{1-c_1-c_2}{z-\gamma}+\frac{1-d_1-d_2}{z-\delta}\right)\frac{\mathrm d w}{\mathrm dz}+\left(\frac{a_1 a_2\, p^\prime(\alpha)}{z-\alpha}+\frac{b_1 b_2\, p^\prime(\beta)}{z-\beta}+\frac{c_1 c_2\, p^\prime(\gamma)}{z-\gamma}+\frac{d_1 d_2\, p^\prime(\delta)}{z-\delta}\right)\frac{w}{p(z)}=0$$
(working off the appearance of the more familiar Riemann DE). Here, $p(z)=(z-\alpha)(z-\beta)(z-\gamma)(z-\delta)$, and we have the additional restriction that $a_1+a_2+b_1+b_2+c_1+c_2+d_1+d_2=2$. (Assuming I am understanding Fuchs's theory correctly.)
We want a Möbius transformation that maps $[\alpha,\beta,\delta]\rightarrow[0,1,\infty]$, which leads to the Möbius transformation $t=\dfrac{z-\alpha}{\beta-\alpha}\dfrac{\beta-\delta}{z-\delta}$. Additionally, we need to make the substitution $w(z)=\left(\dfrac{z-\alpha}{z-\delta}\right)^{a_1} \left(\dfrac{z-\beta}{z-\delta}\right)^{b_1} \left(\dfrac{z-\gamma}{z-\delta}\right)^{c_1}y(t)$.
After much algebra, sweat, and tears, I find that the solution to the original Fuchsian ODE should now look like this:
$$\left(\frac{z-\alpha}{z-\delta}\right)^{a_1} \left(\frac{z-\beta}{z-\delta}\right)^{b_1} \left(\frac{z-\gamma}{z-\delta}\right)^{c_1}\\ \operatorname{\mathit{H\ell}}\left(\dfrac{\gamma-\alpha}{\beta-\alpha}\dfrac{\beta-\delta}{\gamma-\delta},q;a_1+b_1+c_1+d_1,a_1+b_1+c_1+d_2,1+a_1-a_2,1+b_1-b_2;\frac{z-\alpha}{\beta-\alpha}\dfrac{\beta-\delta}{z-\delta}\right)$$
Herein lies my problem: the expression for the accessory parameter $q$ seems to look absolutely horrendous:
$$q=-\frac{(a_1+b_1) (a_2+b_2-1) (\beta -\delta)^2}{(\alpha -\beta) (\gamma -\delta)}+\frac{a_1 a_2 (\gamma -\alpha)+(\beta -\delta) (-a_1 a_2-a_1 b_2+a_1-a_2 b_1+b_1)}{\gamma -\delta}+\frac{(\beta -\delta) (a_1 a_2+a_1 b_2-a_1+a_2 b_1+b_1 b_2-b_1+c_1 c_2-d_1 d_2)}{\alpha -\beta}-a_1 a_2-a_1 c_2+a_1-a_2 c_1-\frac{c_1 c_2 (\gamma -\delta)}{\alpha -\beta}+c_1$$
My already very limited algebraic abilities, having been sufficiently taxed by previous manipulations, was of no avail in finding a simpler/more elegant expression for $q$. Does it really look like this, or might there be a manipulation/substitution I am missing that will result in a tidy-looking expression for $q$? (Once more, I give the reminder that $a_1+a_2+b_1+b_2+c_1+c_2+d_1+d_2=2$; I attempted to use this relation in the initially even more horrendous expression for $q$ that came out of manipulating the Fuchsian ODE, resulting in the only slightly less awkward expression shown above.)
On that note: did I do all the right substitutions, or were there any other substitutions I should have done instead that would have made the algebra less burdensome?
(This and this are definitely related questions, both answered by Alexandre Eremenko, but I do not have access to the references pointed out in the first linked thread, so I am unable to confirm.)
