WriteHeun := proc(F, k, l, m)
	local e0,e1,ei,E,Dx,x,Le,L,v,T,i,j,S1;
	if not assigned(_Envdiffopdomain) then
		error "Should assign value to _Envdiffopdomain"
	fi;
	Dx,x := op(_Envdiffopdomain);

	Le := (x-1)*x*Dx^2+(-1+2*x+e0-e0*x-e1*x)*Dx+1/4*(e0-1+e1+ei)*(-1+e0+e1-ei);
	# has exponents [0,e0] at x=0,  [0,e1]  at x=1,  and exponent-difference ei at x=infinity.

	E := {e0 = 1/l, e1 = 1/k, ei = 1/m};
	L := eval(Le, E);

	# Apply change of variables x --> F:
	L := add(DEtools[mult]( subs(x=F,coeff(L,Dx,i)) , (1/diff(F,x) * Dx)$i),i=0..2);
	L := collect(L/lcoeff(L,Dx), Dx, evala);

	# Now eliminate false singularities with exponent-difference 1, for instance,
	# if the exponents are e, e+1 at x=p then multiply the solutions by (x-p)^(-e)
	# to obtain exponents 0, 1.

	S1 := coeff(L,Dx,1)/2;
	L := DEtools[symmetric_product](L, Dx - coeff(L,Dx,1)/2);
	v := factors(denom(evala(coeff(L, Dx, 0))), indets(f,{RootOf,radical,nonreal}))[2];
	v := add(min(seq(j[1], j=DEtools[gen_exp](L,T,x=RootOf(i[1],x)))) * diff(i[1],x)/i[1], i=v);
	S1 := evala(S1 - v);
	L := evala(Primpart( DEtools[symmetric_product](L, Dx + v) , Dx));
	L := collect(L,Dx,factor);

	sort(L,Dx),  MakeTcoeff1(`DEtools/kovacicsols/e_int`(S1, x)) *  eval(
		hypergeom([1-e0-e1-ei, 1-e0-e1+ei]/2 ,[1-e0],F), E)
end;

# Adjust an expression by a constant factor to make it easier for "series" to process the expression:
MakeTcoeff1 := proc(f)
	local x;
	x := _Envdiffopdomain[2];
	if not has(f,x) then
		f
	elif type(f,`*`) then
		map(procname,f)
	elif type(f,`^`) then
		procname(op(1,f))^op(2,f)
	elif type(f,polynom(anything,x)) then
		collect(f/tcoeff(f,x),x)
	else
		error "unexpected"
	fi
end:


# Example B34 from http://www.math.fsu.edu/~hoeij/Heun/BelyiMaps

_Envdiffopdomain := [Dx,x]:
B34 := 16/27*(128*x-103)*(128*x^7-511*x^6-1232*x^5+3570*x^4+4480*x^3+33208*x^2+56224*x-22972)^3/(8*x+15)/(7*x-32)^2/(19*x^2-40)^7;

# In the notation of the paper we have (k,l,m) with k <= l <= m and the exponent-differences
# of the hypergeometric equation at 1, 0, infinity (in that order) are 1/k, 1/l, 1/m.

WriteHeun(B34, 2, 3, 7);  # B34 has k, l, m = 2, 3, 7


# You can download this file at: www.math.fsu.edu/~hoeij/Heun/StandardForm
read "/home/m1/hoeij/public_html/Heun/StandardForm":

# This tries to move points to 0,1,infinity,t (if that is possible without a field extension) to end up with a Heun equation

_Env_StandardForm := 2;
V := StandardForm(B34);
F := V[1];

WriteHeun(F, 2, 3, 7);

# Lets swap 0,infinity so that the factor x appears in the numerator (that also means swapping l,m)
# That way we can apply "series" to Sol.
F := 1/F;
EQ, Sol := WriteHeun(F, 2, 7, 3);

# Check that Sol is indeed a solution of the Heun equation EQ:
EQ := DEtools[diffop2de](EQ, y(x));
simplify( eval(EQ, y(x) = Sol) );


# An example of Heun( .. ) = .. * hypergeom( .. )
#
coeff(rhs(dsolve(EQ)), _C1) = Sol;

# Check this last equation:
series(lhs(%) - rhs(%), x=0, 20);