//--------------------------Calculating Heun equations from 38 examples of Herfurtner------------------------
LIB "foliation.lib"; LIB "linalg.lib"; LIB "painleve-heun.lib"; ring r=(0,a,b,c,z,rh),x,dp; 
//a,b,c are used for possible exponent parameters, z is the  variable of differential equation and
// rh is the possible extra parameter. 
//Here you can put any example in the Herfurtner list.
number g2=3*(16*z^4-16*z^2+1); number g3=64*z^6-96*z^4+30*z^2+1;
//This command checks factorization of f(x) 
list gzd=newgzweigdrei(g2,g3); g2=gzd[1]; g3=gzd[2];
if (gzd[3]==1) {matrix AAA=gzweigdrei(g2,g3,z,a);}
if (gzd[3]==2) {matrix AAA=gzweigdrei(g2,g3,z,a,b);}
if (gzd[3]==3) {matrix AAA=gzweigdrei(g2,g3,z,a,b,c);}
list lll=lineardiff(AAA,z); number sss=SLForm(lll[1],lll[2],z);
//From this point on, we have to treat 1. families with rational singularities 
//2. families with Galois conjugated singularities 3. 21th family, separately.
list singo=1,0,-1; painheunpar(sss,z,singo); //For Example 1 we run
number g2=3*z*(9*z^3-8); number g3=27*z^6-36*z^3+8;
g2=number(substpar(g2,z,z+1)); g3=number(substpar(g3,z,z+1));
list gzd=newgzweigdrei(g2,g3); g2=gzd[1]; g3=gzd[2];
if (gzd[3]==1) {matrix AAA=gzweigdrei(g2,g3,z,a);}
if (gzd[3]==2) {matrix AAA=gzweigdrei(g2,g3,z,a,b);}
if (gzd[3]==3) {matrix AAA=gzweigdrei(g2,g3,z,a,b,c);}
list lll=lineardiff(AAA,z); number sss=SLForm(lll[1],lll[2],z);
heunparq(sss,z,3,3); //the roots of z^2+3z+3=0 are two Galois conjugated singularities.