Fuchsian systems associated to families of elliptic curves with $4$ singular fiber:

S. Herfurtner 1991 classified families of elliptic curves with exactly four singular fibers. Among them only five families depend on an extra parameter such that the j-function of the family is not constant with respect to this extra parameter. Ch. Doran 2001 explained how to use these five families and obtain algebraic solutions of sixth Painleve equation. Gavrilov and Ben Hamed 2005 explained how all these is related to common algebraic solutions of sixth Painleve equations. Together with Stefan Reiter, we have calculated the Fuchsian system corresponding to to the five families in Herfurtner list. We have also classified the parameters in which the corresponding Painleve equation has an algebraic solution comming from geometry, i.e it comes from the Gauss-Manin connection of families of algebraic varieties. Bellow is the commands used in Singular and the corresponding information for five families. 
LIB "foliation1.34.lib"; LIB "linalg.lib";
ring r=(0,s,a,b,c,z,t_2,t_3),x,dp;
//s is a parameter used inside two procedures in http://w3.impa.br/~hossein/foliation-allversions/ramanujan.txt
//a is the parameter of y=(..)^a and b is the parameter of the deformable family of elliptic curves.
//c is needed when we can write y=(..)^a(...)^c (just for second example in the Herfurtner list). 
// z is the parameter of the family of elliptic curves. For 4 values of z we have singular elliptic curve.
// t_2, t_3 are parameters of an elliptic curve in the Weierstrass form y^2=4x^3-t_2x-t_3

//--------------------The system in variables t_2, t_3
//****We run the two procedures in http://w3.impa.br/~hossein/foliation-allversions/ramanujan.txt *******
//*******************************************************************************************************
//For the case y=(4x^3-t_2x-t_3)^a we run
     poly g=4*x^3-t_2*x-t_3;
     poly gf= a*diff(g,x);
     poly gtf_2= a*diffpar(g,t_2);
     poly gtf_3= a*diffpar(g,t_3);

     matrix gm2[2][2]=gaussmaninmatrix(t_2,g, gf,gtf_2);
     matrix gm3[2][2]=gaussmaninmatrix(t_3,g, gf,gtf_3);
     // to make the fundamental system with det=const.
     matrix ide[2][2]=1,0,0,1;
     gm2=gm2-((1/2)*(1/2-a)*3*(-t_2^2)/(27*t_3^2-t_2^3))*ide;
     gm3=gm3-((1/2)*(1/2-a)*54*t_3/(27*t_3^2-t_2^3))*ide; 
//For the case y=(4x^2-t_2x+t_3)^c(x+t_2/4)^a we run
    poly g=(x+t_2/4)*(4*x^2-t_2*x+t_3);
    poly gf= a*(4*x^2-t_2*x+t_3)+c*(8*x-t_2)*(x+t_2/4);
    poly gtf_2=(c*(x+t_2/4)*(-x)+a*(1/4)*(4*x^2-t_2*x+t_3));
    poly gtf_3= c*(x+t_2/4);

    matrix gm2[2][2]=gaussmaninmatrix(t_2,g, gf,gtf_2);
    matrix gm3[2][2]=gaussmaninmatrix(t_3,g, gf,gtf_3);
    // to make the fundamental system with det=const.
    matrix ide[2][2]=1,0,0,1;
    gm2=gm2-1/2*((1-c-a)*2*t_2/(t_2^2+2*t_3)+(1/2-c)*2*t_2/(t_2^2-16*t_3) )*ide;
    gm3=gm3-1/2*((1-c-a)*2/(t_2^2+2*t_3)+(1/2-c)*(-16)/(t_2^2-16*t_3) )*ide;
// For the rest we only need the matrices gm2, gm3.--------------------------------------------------------
//*********************************************************************************************************
//For the case y=(4x^3-t_2x-t_3)^a we have
matrix gm2[2][2]=(-1/4*t_2^2)/(t_2^3-27*t_3^2),(27*a*t_3-18*t_3)/(t_2^3-27*t_3^2), (9/4*a*t_2*t_3-3/4*t_2*t_3)/(t_2^3-27*t_3^2), (1/4*t_2^2)/(t_2^3-27*t_3^2);
matrix gm3[2][2]=(9/2*t_3)/(t_2^3-27*t_3^2), (-18*a*t_2+12*t_2)/(t_2^3-27*t_3^2), (-3/2*a*t_2^2+1/2*t_2^2)/(t_2^3-27*t_3^2), (-9/2*t_3)/(t_2^3-27*t_3^2);
//For the case y=(4x^2-t_2x+t_3)^c(x+t_2/4)^a we have
matrix gm2[2][2]=(6*a*t_2*t_3-6*c*t_2*t_3-1/2*t_2^3+5*t_2*t_3)/(t_2^4-14*t_2^2*t_3-32*t_3^2), (-48*a*t_3-96*c*t_3+96*t_3)/(t_2^4-14*t_2^2*t_3-32*t_3^2), (12*a*t_3^2-3*c*t_2^2*t_3+t_2^2*t_3-4*t_3^2)/(t_2^4-14*t_2^2*t_3-32*t_3^2), (-6*a*t_2*t_3+6*c*t_2*t_3+1/2*t_2^3-5*t_2*t_3)/(t_2^4-14*t_2^2*t_3-32*t_3^2);
matrix gm3[2][2]=(-3*a*t_2^2+3*c*t_2^2+t_2^2+8*t_3)/(t_2^4-14*t_2^2*t_3-32*t_3^2), (24*a*t_2+48*c*t_2-48*t_2)/(t_2^4-14*t_2^2*t_3-32*t_3^2), (-12*a*t_2*t_3+3*c*t_2^3-t_2^3+4*t_2*t_3)/(2*t_2^4-28*t_2^2*t_3-64*t_3^2), (3*a*t_2^2-3*c*t_2^2-t_2^2-8*t_3)/(t_2^4-14*t_2^2*t_3-32*t_3^2);
//***********************************************************************************************************
//---------------------------------------Herfurtner list with parameter-----------------------------
//***********************************************************************************************************
//We have to substitute b for somthing so that all the singularities becomes polynomial expressions in b
// 5 list with only a parameter
//list 1, I1-I1-II-IV*
poly g2=3*(z-1)*(z-b^2)^3;
poly g3=(z-1)*(z-b^2)^4*(z+b);

//list 2, I1-I1-I2-I2*
poly bb=3/4*(b+1/b)+1/2;
poly g2=12*z^2*(z^2+bb*z+1);
poly g3=4*z^3*(2*z^3+3*bb*z^2+3*bb*z+2);

//list 3  I1-I1-I1-I3*
poly bb=2/3*(b+1/b)-1/3;
poly g2=12*z^2*(z^2+2*bb*z+1);
poly g3=4*z^3*(2*z^3+3*(bb^2+1)*z^2+6*bb*z+2);
poly sinset=12*x^2+3*(3*bb^2+6*bb-1)*x+4*(bb+2);

//list 4  I1-I1-I1-III*
poly bb=(2/3)*(((b^2)/3-1)/((b^2)/3+1))+1/3;
poly g2=3*z^3*(z+bb);
poly g3=z^5*(z+1);
poly sinset=(3*bb-2)*x^2+(3*bb^2-1)*x+bb^3;

//list 5 I1-I1-I2-IV*
poly bb=1/4*(b+2/b);
poly g2=3*z^3*(z+2*bb);
poly g3=z^4*(z^2+3*bb*z+1);
poly sinset=(3*bb^2-2)*x^2+2*bb*(4*bb^2-3)*x-1;
//-----------------------The list 2 I1-I1-I2-I2* with parameters a and c
poly bb=3/4*(b+1/b)+1/2;
//poly g2=12*z^2*(z^2+bb*z+1);
//poly g3=4*z^3*(2*z^3+3*bb*z^2+3*bb*z+2);
//In order to identify g2 and g3 we had to run  factorize(4*x^3-g2*x-g3);
poly g2=4*(z^2+z); 
poly g3=(-9*b^2*z^3-8*b*z^4+2*b*z^3-8*b*z^2-9*z^3)/(b);
//this g2 and g3 represents t2 and t3
//--------------------------------------The system after pullback in z variable
matrix komak2=substpar(gm2,t_2,g2); komak2=substpar(komak2, t_3,g3); 
matrix komak3=substpar(gm3,t_2,g2); komak3=substpar(komak3, t_3,g3);
matrix sm=komak2*diffpar(g2,z)+komak3*diffpar(g3,z);
//-------------------------------------To write in the SL form of the first coordinate
number hilfe=diffpar(number(sm[1,2]),z)/number(sm[1,2]);
//The Picard Fuchs equation y''+p1y'+p2y=0 of the first coordinate of the system  
number p1=-hilfe-number(sm[1,1])-number(sm[2,2]);
number p2=number(det(sm))-diffpar(number(sm[1,1]),z)+number(sm[1,1])*hilfe;
//SL form y''=slform*y
number slform=-(p2-(1/2)*diffpar(p1,z)-(1/4)*p1^2);
//-------------------------Here we have to find the apparant singularity. We have used the commands
                            //poly  poles=denominator(slform);
                            //poles=substpar(poles, z,x);
                            //factorize(poles);
                            //factorize(sinset);
//************************************************************************
// 5 list with only a parameter
//The second singularity sing2 must be zero
//list 1
number appa=b;
number sing1=b^2;
number sing2=0;
number sing3=1;

//list 2
number appa=1;
number sing1=-b;
number sing2=0;
number sing3=-1/b;

//list 3
number appa=-(-2*b^2-5*b-2)/(3*b);
number sing1=-(2*b+1)/(3*b^2);
number sing2=0;
number sing3=-(b^2+2*b)/3;

//list 4
number appa=-(b^2-1)/(b^2-9);
number sing1= -(b^3+3*b^2+3*b+1)/(b^3+3*b^2+3*b+9); 
number sing2=0;
number sing3=-(b^3-3*b^2+3*b-1)/(b^3-3*b^2+3*b-9);

//list 5
number appa=-(-4*b^3-8*b)/(3*b^4-20*b^2+12);
number sing1=8/(b^3-6*b);
number sing2=0;
number sing3=-(2*b^3)/(3*b^2-2);
//-------------------------
//*****************************************************************************
//--The SL form written in Yoshida's book p. 173. Note that 1 may not be a singular point
number aeins, azwei, adrei, avier, afunf, nue, ele;  
aeins=number(substpar(slform*(z-sing1)^2, z, sing1));
azwei=number(substpar(slform*(z-sing2)^2, z, sing2));
adrei=number(substpar(slform*(z-sing3)^2, z, sing3));
afunf=number(substpar(slform*(z-appa)^2, z, appa));
if (afunf<>3/4){"There are some problems";}

slform=slform-aeins/(z-sing1)^2-azwei/(z-sing2)^2-adrei/(z-sing3)^2-afunf/(z-appa)^2;


nue=-number(substpar(slform*(z-appa), z, appa))*sing3;
ele=number(substpar(slform*(z-sing1), z, sing1))*sing3;
slform=slform+appa*(appa-sing3)*nue/((z-sing2)*(z-sing3)*(z-appa)*sing3)-sing1*(sing1-sing3)*ele/((z-sing2)*(z-sing3)*(z-sing1)*sing3);

avier=slform*(z-sing2)*(z-sing3);
//Calculating the exponents, see Yoshida p. 174
number the1h, the2h, the3h, the4h;
number alpha, beta, gamma, delta;

the1h=4*aeins+1;
the2h=4*azwei+1;
the3h=4*adrei+1;
the4h=(the1h+the2h+the3h-1)+4*avier+2;
alpha=(1/2)*the4h;
beta=(1/2)*the2h;
gamma=(1/2)*the3h;
delta=(1/2)*the1h;
//**************************************************************************************************
//---------------------The exponents of the examples
//list 1
number the1=1/3;
number the2=a-1/2;
number the3=1/3;
number the4=(a-1/2);

//list 2
number the1=a-1/2;
number the2=2*a-1;
number the3=a-1/2;
number the4=2*a-1;

//list 3
number the1=a-1/2;
number the2=3*(a-1/2);
number the3=a-1/2;
number the4=a-1/2;

//list 4
number the1=a-1/2;
number the2=1/2;
number the3=a-1/2;
number the4=a-1/2;

//list 5
number the1=a-1/2;
number the2=1/3;
number the3=a-1/2;
number the4=2*(a-1/2);

//list 2 with parameters a and c
number the1=c-1/2;
number the2=a+c-1;
number the3=c-1/2;
number the4=a+c-1;

//**************************************************************************************************
//------Finding the apparant singularity and movable singularity Yoshida p.174
number lambda=appa/sing3;
number tparam=sing1/sing3;


//a mixture of Yoshida p. 174 (4.2.5) (4.2.6)
number Mu=nue-(1/2)*((1-the1)/(lambda-tparam)+(1-the2)/lambda+(1-the3)/(lambda-1));
//Yoshida p. 208
number Meins=(lambda-tparam)/(tparam*(tparam-1));
number Mzwei=lambda/tparam;
number Mdrei=(lambda-1)/(1-tparam);
number Meinseins=lambda*(lambda-1);
//Yoshida p.
number Alph=(-1/2)*(the1+the2+the3+the4-1);
number Weins=Meinseins*(Mu+Alph/lambda);
number Wzwei=(1/Mzwei)*((tparam-1)*Meins*Weins-Alph);
number Wdrei=(-1/Mdrei)*tparam*Meins*Weins;
number WW=(1/(the4-1))*(Weins*(Meins*Weins-the1)+Wzwei*(Mzwei*Wzwei-the2)+Wdrei*(Mdrei*Wdrei-the3));


list Mlist=Meins,Mzwei, Mdrei;
list Wlist=Weins,Wzwei, Wdrei;
list Tlist=the1,the2,the3;
//--------Garnier to Schlesinger Yoshida p. 208
matrix Qeins[2][2];
matrix Qzwei[2][2];
matrix Qdrei[2][2];
list QQ=Qeins, Qzwei, Qdrei;
for (int i=1; i<=3; i=i+1)
    {
    QQ[i][1,1]=Mlist[i]*(Wlist[i]-WW);
    QQ[i][1,2]=-Mlist[i];
    QQ[i][2,1]=-(Wlist[i]-WW)*(Mlist[i]*(WW-Wlist[i])+Tlist[i]);
    QQ[i][2,2]=Tlist[i]-Mlist[i]*(Wlist[i]-WW);
    }

matrix fsm=(1/(z-tparam))*QQ[1]+(1/z)*QQ[2]+(1/(z-1))*QQ[3];

//--------------------------
matrix coord[1][2]=0,1;
poly che=sysdif(fsm, coord,z)[1,3];//We can use directly the formula
che=che/(z*(z-1)*(z-tparam));
number nlambda=-number(substpar(che,z,0)/diffpar(che,z));
//Okamato transformation obtained by permuting the coordinates of a solution
the4=2-the4;
alpha=(1/2)*(the4)^2;
//-------------------------------------------------
//******If you want to check the solution lambda, tparam run this; number nlambda=lambda; ********************
//------------To check that the obtained algebraic solutions are really solutions.----------------------
number nlambdap=diffpar(nlambda,b); number tparamp=diffpar(tparam,b); 
number nlambdatt=((diffpar(tparam,b))*diffpar(nlambdap,b)-nlambdap*diffpar(tparamp,b))/(tparamp)^3;
number nlambdat=nlambdap/(tparamp);
number painleve1=(1/2)*(1/nlambda+1/(nlambda-1)+1/(nlambda-tparam))*nlambdat^2;
number painleve2=-(1/tparam+1/(tparam-1)+1/(nlambda-tparam))*nlambdat;
number painleve31=nlambda*(nlambda-1)*(nlambda-tparam)/(tparam^2*(tparam-1)^2);
number painleve32=alpha-beta*tparam*(1/(nlambda^2))+gamma*(tparam-1)*(1/((nlambda-1)^2))+(1/2-delta)*tparam*(tparam-1)*
(1/((nlambda-tparam)^2));


number painleve=nlambdatt-painleve1-painleve2-painleve31*painleve32;