**Appendix A**

The following Mathematica package implements the algorithm producing the coefficients of the method for double precision. In the input we provide the free coefficients *c*2, *c*4, *c*5, *c*6, *c*<sup>7</sup> and ˆ *b*9. In the output we get four vectors, namely *b*, ˆ *b*, *c* and the matrix *A*.

```
NEW65[cc2_, cc4_, cc5_, cc6_, cc7_, bbb9_] :=
Module[{b, a, c, bb, vandh, vandl, ac, ac2, cc, ii, ba, cond, soh, b1, b4, b5, b6,
b7, b8, bb1, bb4, bb5, bb6, bb7, bb8, bb9, c2, c3, c4, c5, c6, c7, a21, a31,
a32, a41, a43, a51, a53, a54, a61, a63, a64, a65, a71, a73, a74, a75, a76,
a81, a83, a84, a85, a86, a87, so3, so5, sol, so6, so7, so8},
c2 = Rationalize[cc2, 10^-16]; c4 = Rationalize[cc4, 10^-16];
c5 = Rationalize[cc5, 10^-16]; c6 = Rationalize[cc6, 10^-16];
c7 = Rationalize[cc7, 10^-16]; bb9 = Rationalize[bbb9, 10^-16];
b={b1,0,0,b4,b5,b6,b7,b8,0};
a={{0,0,0,0,0,0,0,0,0},
{a21,0,0,0,0,0,0,0,0},
```

```
{a31,a32,0,0,0,0,0,0,0},
{a41,0,a43,0,0,0,0,0,0},
{a51,0,a53,a54,0,0,0,0,0},
{a61,0,a63,a64,a65,0,0,0,0},
{a71,0,a73,a74,a75,a76,0,0,0},
{a81,0,a83,a84,a85,a86,a87,0,0},
{b1,0,0,b4,b5,b6,b7,b8,0}};
c={0,c2,c3,c4,c5,c6,c7,1,1};
bb={bb1,0,0,bb4,bb5,bb6,bb7,bb8,bb9};
vandh={b.c==1/2,b.c^2==1/3,b.c^3==1/4,b.c^4==1/5,b.c^5==1/6};
vandl={bb.c==1/2,bb.c^2==1/3,bb.c^3==1/4,bb.c^4==1/5};
ac=a.c-c^2/2;
ac2=a.c^2-c^3/3; cc=DiagonalMatrix[c]; ii=IdentityMatrix[9];
ba=b.(a+cc-1*ii);
cond=
{b.(cc-1*ii).a.(cc-c4*ii).(cc-c5*ii).c-Integrate[(x-1)*Integrate[(x-c5)*(x-c4)*x,
{x,0,x}],{x,0,1}],
b.(cc-ii).a.(cc-c4*ii).(cc-c5*ii).c-Integrate[(x-1)*Integrate[(x-c4)*(x-c5)*x,
{x,0,x}],{x,0,1}] /.{a43->0,a53->0,a63->0,a73->0},
bb.a.(cc-c5*ii).(cc-c4*ii).c-Integrate[Integrate[(x-c5)*(x-c4)*x,{x,0,x}],
{x,0,1}]};
(* start procedure *)
soh=Solve[vandh,{b4,b5,b6,b7,b8}];
b4=Together[soh[[1,1,2]]];
b5=Together[soh[[1,2,2]]];
b6=Together[soh[[1,3,2]]];
b7=Together[soh[[1,4,2]]];
b8=Together[soh[[1,5,2]]]; (* OK b *)
sol=Solve[Join[{(bb.a)[[3]]==0},vandl],{bb4,bb5,bb6,bb7,bb8}];
bb4=Together[sol[[1,1,2]]];
bb5=Together[sol[[1,2,2]]];
bb6=Together[sol[[1,3,2]]];
bb7=Together[sol[[1,4,2]]];
bb8=Together[sol[[1,5,2]]]; (* OK bb *)
c3=2/3*c4; a43=c4^2/(2*c3); a32=c3^2/(2*c2);
so5=Solve[{ac[[5]]==0,ac2[[5]]==0},{a53,a54}];
a53=Together[so5[[1,1,2]]];a54=Together[so5[[1,2,2]]]; (* after upper left *)
a87=Together[Solve[{ba[[7]]==0},{a87}][[1,1,2]]];
a76=Together[Solve[{cond[[2]]==0},{a76}][[1,1,2]]];
a86=Together[Solve[{ba[[6]]==0},{a86}][[1,1,2]]];
ac=Together[ac]; ac2=Together[ac2];
ba=Together[ba]; cond=Together[cond]; (* OK down right *)
so3=Solve[{ba[[3]]==0,cond[[1]]==0,cond[[3]]==0},{a63,a73,a83}];
a63=Together[so3[[1,1,2]]];
a73=Together[so3[[1,2,2]]];
a83=Together[so3[[1,3,2]]]; (* OK 3nd column *)
so6=Solve[{ac[[6]]==0,ac2[[6]]==0},{a64,a65}];
a64=Together[so6[[1,1,2]]];a65=Together[so6[[1,2,2]]];
so7=Solve[{ac[[7]]==0,ac2[[7]]==0},{a74,a75}];
a74=Together[so7[[1,1,2]]];a75=Together[so7[[1,2,2]]];
so8=Solve[{ac[[8]]==0,ac2[[8]]==0},{a84,a85}];
a84=Together[so8[[1,1,2]]];a85=Together[so8[[1,2,2]]]; (* OK ai4, ai5 *)
b1=1-b4-b5-b6-b7-b8;
bb1=1-bb4-bb5-bb6-bb7-bb8;
a21=c2;
a31=c3-a32;
a41=c4-a43;
a51=c5-a53-a54;
a61=c6-a63-a64-a65;
```
a71=c7-a73-a74-a75-a76; a81=1-a83-a84-a85-a86-a87; Return[SetAccuracy[{b, bb, c, a}, 16] // Chop]]
