BeginPackage[ "csa`"];
Clear[ "csa`*" ]

TEST::usage = "  TEST[x1,x2,x3,x4]  of a RK5(4) with b.(a+c-i)=0 "

Begin["`Private`"];
Clear[ "csa`Private`*" ];

TEST[cc2_?NumericQ, cc3_?NumericQ, cc4_?NumericQ, bc2_?NumericQ] :=
  Module[{a32, a42, a43, a52, a53, a54, a62, a63, a64, a65, bb1, bb2, bb3, bb4, bb5, bb6, bb7, b1, b2, b3, b4, b5, b6, 
          c2, c3, c4, c5, so, so1, a, b, bb, c, ba, equ, eequ, temp, j1, tec6}, 

c2 = Rationalize[cc2, 10^-16]; c3 = Rationalize[cc3, 10^-16]; c4 = Rationalize[cc4, 10^-16];b2 = Rationalize[bc2, 10^-16];

equ = {-(1/120) + b.a.a.a.c, -(1/60) + b.a.a.c^2, -(1/40) + b.a.(c*a.c),
       -(1/20) + b.a.c^3, -(1/30) + b.(c*a.a.c), -(1/15) +  b.(c*a.c^2),
       -(1/20) + b.(a.c)^2, -(1/10) + b.(c^2*a.c), -(1/5) +  b.c^4,
       -(1/24) + b.a.a.c, -(1/12) + b.a.c^2, -(1/8) + b.(c*a.c),
       -(1/4) + b.c^3, -(1/6) + b.a.c, -(1/3) + b.c^2, -(1/2) + b.c,
       -1 + b.{1, 1, 1, 1, 1, 1, 1}};

eequ = {-(1/24) + bb.a.a.c, -(1/12) + bb.a.c^2, -(1/8) +  bb.(c*a.c),
        -(1/4) + bb.c^3, -(1/6) + bb.a.c, -(1/3) + bb.c^2, -(1/2) + bb.c,
        -1 + bb.{1, 1, 1, 1, 1, 1, 1}};

c = {0, c2, c3, c4, c5, 1, 1};

b = {b1, b2, b3, b4, b5, b6,  0};

bb = {bb1, bb2, bb3, bb4, bb5, bb6, 1/40};

a = {{0, 0, 0, 0, 0, 0, 0},
     {c2, 0, 0, 0, 0, 0, 0},
     {c3-a32, a32, 0, 0, 0, 0, 0},
     {c4 - a42 - a43, a42, a43, 0,  0, 0, 0},
     {c5 - a52 - a53 - a54, a52, a53, a54, 0, 0,  0},
     {1 - a62 - a63 - a64 - a65, a62, a63, a64, a65, 0, 0}, b};

so = Solve[{equ[[17]] == 0, equ[[16]] == 0, equ[[15]] == 0, equ[[13]] == 0, equ[[9]] == 0}, {b1, b3, b4, b5, b6}];

{b1, b3, b4, b5, b6} = Simplify[so[[1, 1 ;; 5, 2]]];
so = Solve[{eequ[[4]] == 0, eequ[[6]] == 0, eequ[[7]] == 0, eequ[[8]] == 0}, {bb1, bb3, bb4, bb5}];

{bb1, bb3, bb4, bb5} = Simplify[so[[1, 1 ;; 4, 2]]];

ba = Simplify[b.(a + DiagonalMatrix[c] - IdentityMatrix[7])];

so = Solve[{ba[[2 ;; 5]] == {0, 0, 0, 0}}, {a62, a63, a64, a65}];

{a62, a63, a64, a65} = Simplify[so[[1, 1 ;; 4, 2]]];

so = Solve[{equ[[2]] == 0, equ[[1]] == 0, equ[[3]] == 0, equ[[8]] == 0}, {a43, a52, a53, a54}, Sort -> False];

{a43, a52, a53, a54} = Simplify[so[[1, 1 ;; 4, 2]]];

{a62, a63, a64, a65} = Simplify[{a62, a63, a64, a65}];

so = Solve[{eequ[[1]] == 0, eequ[[2]] == 0}, {bb2, bb6}, Sort->False];

{bb2, bb6} = Simplify[so[[1, 1 ;; 2, 2]]];

so=Solve[{eequ[[5]]==0},{a42}];a42=Simplify[so[[1,1,2]]];

{a43, a52, a53, a54, a62, a63, a64, bb1, bb3, bb4, bb5} =
                    Simplify[{a43, a52, a53, a54, a62, a63, a64, bb1, bb3, bb4, bb5}];

so= Solve[{eequ[[3]] == 0}, {c5}]; 
c5 = Simplify[so[[1, 1, 2]]];

{a52, a53, a54, a43, bb5, bb6, a42, b1, b2, b3, b4, b5, b6} = Simplify[{a52, a53, a54, a43, bb5, bb6, a42, b1, b2, b3, b4, b5, b6}];

{bb1, bb2, bb3, bb4, bb5, bb6} = Simplify[{bb1, bb2, bb3, bb4, bb5 ,bb6}];

so = Union[Solve[{Numerator[Simplify[equ[[7]]]] == 0}, {a32}]];

temp = N[Select[so, Im[#[[1, 2]]] == 0 &],25];

Return[Table[{a, b, bb, c} /. temp[[j1]], {j1, 1, Length[temp]}]]
];

End[ ];
EndPackage[ ];

(* --------------------------------- *)
(* ------ The proposed pair -------- *)
TEST[0.161,0.327,0.9,0.01][[6]]

a={{0, 0, 0, 0, 0, 0, 0}, 
   {0.161, 0, 0, 0, 0, 0, 0}, 
   {-0.00848065549235698854, 0.335480655492356989, 0, 0, 0, 0, 0}, 
   {2.89715305710549343, -6.35944848997507484, 4.36229543286958141, 0, 0, 0, 0}, 
   {5.325864828439257, -11.748883564062828, 7.49553934288983621, -0.09249506636175525, 0, 0, 0}, 
   {5.8614554429464200, -12.9209693178471093, 8.1593678985761586, -0.071584973281400997, -0.0282690503940683829, 0, 0}, 
   {0.0964607668180652295, 0.01, 0.479889650414499575, 1.37900857410374189, -3.2900695154360807, 2.32471052409977398, 0}};
b= {0.0964607668180652295, 0.01, 0.479889650414499575, 1.37900857410374189, -3.2900695154360807, 2.32471052409977398, 0};
bb={0.0935237485818927066,0.00865288314156636761,0.492893099131431868,1.14023541226785810,-2.3291801924393646,1.56887504931661552,0.025};
c= {0, 0.161, 0.327, 0.9, 0.980025540904509686, 1, 1};

(* weights of 4th order interpolant *)
bs={t-2.7637061972748258687 t^2+2.9132554618219126554 t^3-1.0530884977290215572 t^4,
0.13170000000000000000 t^2-0.22340000000000000000 t^3+0.10170000000000000000 t^4,
3.9302962368947515174 t^2-5.9410338721315047357 t^3+2.4906272856512527931 t^4,
-12.411077166933677039 t^2+30.338188630282321650 t^3-16.548102889244902718 t^4,
37.509313416511039179 t^2-88.178904894766401078 t^3+47.379521962819281219 t^4,
-27.896526289197287789 t^2+65.091894674793671508 t^3-34.870657861496609736 t^4,
1.5 t^2-4 t^3+2.5 t^4};