#
# 2013-05-04 Nine stage Runge-Kutta methods of order 7
# (c) 2013 Sergey Khashin
# e-mail: khash2@gmail.com
# www:    http://math.ivanovo.ac.ru/dalgebra/Khashin/index.html
#
#   Free variables c3,c4,c5,a52,c7,c8, BUT ! 
#   for fixex c3,c4,c5,a52 variables c7,c8 saticfys
#   quadratic relation  RK79_c7c8(c3,c4,c5,a52)
#   For rational c3,c4,c5,a52 this relations have raional points
#   which one can find with function ratPoints
#
#
#
# ---- rk_79 ---------------------------------------
# RK79_c6( c3,c4,c5,a52)            c6 as rational function of (c3,c4,c5,a52)
# RK79_tn1 (c3,c4,c5,a52)           service function
# RK79_td1 (c3,c4,c5,a52)           service function
# RK79_L1  (c3,c4,c5,a52)           service function
# RK79_L2  (c3,c4,c5,a52)           service function
# RK79_L4  (c3,c4,c5,a52)           service function
# RK79_b9  (c3,c4,c5,a52,c7,c8)     b9 as rational function of (c3,c4,c5,a52,c7,c8)
# RK79_car1(c3,c4,c5,a52,c7,c8)     -> [b3,b4,b5,b6,b7,b8,b9] for given (c3,c4,c5,a52,c7,c8)
# RK79_c7c8(c3,c4,c5,a52)           biquadratic equation(c7,c8) for given (c3,c4,c5,a52)
# nextX(f,P)                        point P(x,y) -> point(x1,y)
# nextY(f,P)                        point P(x,y) -> point(x,y1)
# point3(f,P1,P2)                   points P1,P2->P3
# ratPoints(f,N::integer)           about 4N rational point on the curve f(x,y)=0
# RK_79a(   c3,c4,c5,a52,c7,c8)     main function, return matrix of 9-stage RK method of order 7
# RK_79b(c2,c3,c4,c5,a52,c7,c8)     main function, return matrix of 9-stage RK method of order 7
#
# RK_79c(   c3,c4,c5,a52,c7,ind)    main function, return matrix of 9-stage RK method of order 7, ind=1|2
# RK_79d(c2,c3,c4,c5,a52,c7,ind)    main function, return matrix of 9-stage RK method of order 7, ind=1|2


# ---------------------- c6 ---------------------------------------------------------------
RK79_c6 := (c3,c4,c5,a52)->
   (112*a52^2*c4*c3^6-432*c4^2*c5^5*c3+252*c4^2*c5^6*c3-36*c4*c5^5*c3-36*c4^2*c5^4*c3
   +1080*c5^4*c3^2*c4^2-378*c5^5*c3^2*c4^2+216*c5^5*c3^2*c4+18*c5^4*c3^2*c4-54*c5^3*c3^2*c4^2
   -126*c5^6*c3^2*c4+27*c5^3*c3^3*c4-378*c5^4*c3^3*c4^2-324*c5^4*c3^3*c4+189*c5^5*c3^3*c4
   -648*c5^3*c3^3*c4^2-12*c5^3*c3^3*a52+54*c5^2*c3^3*c4^2+16*c3^4*a52^2*c4+567*c3^4*c4^2*c5^3
   -12*c3^4*a52*c5^2-192*a52^2*c3^5*c4+112*a52^2*c3^5*c4^2+36*c5^5*c4^2+9*c5^5*c3^2+16*a52^2*c3^5
   +24*c5^3*c3^2*a52*c4+24*c5^2*c3^2*a52*c4^2+12*c5^2*c3^3*a52*c4-288*c5^3*c3^3*a52*c4
   -288*c5^2*c3^3*a52*c4^2+336*c5^3*c3^3*a52*c4^2-24*c5*c3^3*a52*c4^2-252*c3^4*a52*c5^2*c4^2
   -60*c3^4*a52*c5*c4+288*c3^4*a52*c5*c4^2+84*c3^4*a52*c5^3*c4+576*c3^4*a52*c5^2*c4-252*a52*c3^5*c4*c5^2
   -252*a52*c3^5*c4^2*c5)/(108*c5^4*c4^2-1120*a52^2*c4*c3^6+864*c4^2*c5^5*c3-2520*c4^2*c5^6*c3
   +432*c4*c5^5*c3+432*c4^2*c5^4*c3-3672*c5^4*c3^2*c4^2+3780*c5^5*c3^2*c4^2-432*c5^5*c3^2*c4
   -216*c5^4*c3^2*c4+648*c5^3*c3^2*c4^2+1260*c5^6*c3^2*c4-324*c5^3*c3^3*c4+3780*c5^4*c3^3*c4^2
   +1404*c5^4*c3^3*c4-1890*c5^5*c3^3*c4+1296*c5^3*c3^3*c4^2+144*c5^3*c3^3*a52-648*c5^2*c3^3*c4^2
   -192*c3^4*a52^2*c4-5670*c3^4*c4^2*c5^3+144*c3^4*a52*c5^2+832*a52^2*c3^5*c4-1120*a52^2*c3^5*c4^2
   -504*c4*c5^6*c3-108*c4*c5^4*c3-216*c4^2*c5^3*c3+108*c5^3*c3^2*c4+3780*c5^6*c3^2*c4^2
   +108*c5^2*c3^2*c4^2-11340*c5^5*c3^3*c4^2-72*c5^2*c3^3*a52+1134*c3^4*c4^2*c5^2+8505*c3^4*c5^4*c4^2
   -336*c3^4*a52*c5^3+224*c3^4*a52^2*c4^2+224*a52^2*c3^6+1680*a52^2*c3^6*c4^2-432*c5^5*c4^2
   -108*c5^5*c3^2+504*c5^6*c4^2+27*c5^4*c3^2+126*c5^6*c3^2+48*c3^4*a52^2-192*a52^2*c3^5
   -288*c5^3*c3^2*a52*c4-288*c5^2*c3^2*a52*c4^2-144*c5^2*c3^3*a52*c4+912*c5^3*c3^3*a52*c4
   +576*c5^2*c3^3*a52*c4^2-3360*c5^3*c3^3*a52*c4^2+288*c5*c3^3*a52*c4^2+2520*c3^4*a52*c5^2*c4^2
   +720*c3^4*a52*c5*c4-1584*c3^4*a52*c5*c4^2-840*c3^4*a52*c5^3*c4-1152*c3^4*a52*c5^2*c4
   +2520*a52*c3^5*c4*c5^2+2520*a52*c3^5*c4^2*c5+672*c5^3*c3^2*a52*c4^2+144*c5^2*c3^2*a52*c4
   -144*c5*c3^3*a52*c4+5040*c3^4*a52*c5^3*c4^2-1008*a52*c3^5*c4*c5-7560*a52*c3^5*c4^2*c5^2):

# ---------------------- b9 ---------------------------------------------------------------
RK79_tn1 := (c3,c4,c5,a52)->
     -196*a52^2*(136*c4-285*c4^2+180*c4^3-20)*(-1+c5)*c3^7
   +7*a52*(7938*c4*c5^3-6532*c4*a52-9828*c4*c5^2+2520*c4*c5+6532*c4*c5*a52-29295*c4^3*c5^2
      +13088*c4^2*a52-9828*c4^2*c5+7980*c4^3*c5*a52+7938*c4^3*c5+22680*c4^3*c5^3-7980*a52*c4^3
      +1008*a52-1008*c5*a52+37170*c4^2*c5^2-29295*c4^2*c5^3-13088*c4^2*c5*a52)*c3^6
   -1/4*(-436464*c4*c5^2*a52+293412*c4*c5^3*a52+62328*c4*c5^4*a52+1446564*c4^2*c5^2*a52
         -643104*c4^2*c5^3*a52-485100*c4^2*c5^4*a52+341712*c4^3*c5*a52-1083600*c4^3*c5^2*a52
         +396900*c4^3*c5^3*a52+423360*c4^3*c5^4*a52+116592*c4*c5*a52-433608*c4^2*c5*a52
         +79380*c4^2*c5^2-539784*c4^2*c5^3-79380*c4^3*c5^2+539784*c4^3*c5^3+106624*c4^3*c5*a52^2
         -16000*c5*a52^2-96960*c4*a52^2-106624*c4^3*a52^2+182896*c4^2*a52^2+96960*c4*c5*a52^2
         -182896*c4^2*c5*a52^2+10416*c5^2*a52-32928*c5^3*a52+23520*c5^4*a52+16000*a52^2
         +1131165*c4^2*c5^4-714420*c4^2*c5^5-1131165*c4^3*c5^4+714420*c4^3*c5^5)*c3^5
   +1/4*(-146097*c4*c5^4+281232*c4*c5^5-166698*c4*c5^6-176784*c4*c5^2*a52+9696*c4*c5^3*a52
         +134316*c4*c5^4*a52+429984*c4^2*c5^2*a52+486948*c4^2*c5^3*a52-799008*c4^2*c5^4*a52
         +153384*c4^3*c5*a52-276024*c4^3*c5^2*a52-490896*c4^3*c5^3*a52+670320*c4^3*c5^4*a52
         +57840*c4*c5*a52-201360*c4^2*c5*a52+23436*c4*c5^3+119448*c4^2*c5^2-695142*c4^2*c5^3
         -121716*c4^3*c5^2+691362*c4^3*c5^3+15680*c4^3*c5*a52^2-2880*c5*a52^2-16000*c4*a52^2
         -15680*c4^3*a52^2+28224*c4^2*a52^2+16000*c4*c5*a52^2-28224*c4^2*c5*a52^2+14160*c5^2*a52
         -45504*c5^3*a52+32928*c5^4*a52+2880*a52^2+982989*c4^2*c5^4+302400*c4^2*c5^5
         -813645*c4^2*c5^6-894159*c4^3*c5^4-515970*c4^3*c5^5+952560*c4^3*c5^6)*c3^4
   -3/4*c5*(77652*c4*c5^4-3465*c4*c5^5-35574*c4*c5^6-34368*c4*c5^2*a52+32312*c4*c5^3*a52
         +160352*c4^2*c5^2*a52-136920*c4^2*c5^3*a52+21336*c4^3*c5*a52-129080*c4^3*c5^2*a52
         +106624*c4^3*c5^3*a52+1840*c4*c5*a52-22672*c4^2*c5*a52+1440*c5*a52+2880*c4*a52
         -9440*c4^2*a52+540*c5^3+16920*c4^2*c5-17208*c4^3*c5+9540*c4*c5^2-52728*c4*c5^3
         -54180*c4^2*c5^2-120333*c4^2*c5^3+46386*c4^3*c5^2+171486*c4^3*c5^3-4720*c5^2*a52
         +3472*c5^3*a52-3000*c5^4+5292*c5^5-2940*c5^6+6944*a52*c4^3+490686*c4^2*c5^4
         -290304*c4^2*c5^5-74970*c4^2*c5^6-573930*c4^3*c5^4+304290*c4^3*c5^5+105840*c4^3*c5^6)*c3^3
   +3/2*c5^2*(31059*c4*c5^4-21462*c4*c5^5+3472*c4*c5^2*a52-10976*c4^2*c5^2*a52-10976*c4^3*c5*a52
         +7840*c4^3*c5^2*a52-4720*c4*c5*a52+15168*c4^2*c5*a52+1440*c4*a52-4720*c4^2*a52+270*c5^2
         -1500*c5^3+1080*c4*c5+7380*c4^2*c5-8604*c4^3*c5-2730*c4*c5^2-9096*c4*c5^3
         -64128*c4^2*c5^2+101250*c4^2*c5^3+67785*c4^3*c5^2-93954*c4^3*c5^3-1080*c4^3+1080*c4^2
         +2646*c5^4-1470*c5^5+3472*a52*c4^3+8946*c4^2*c5^4-62328*c4^2*c5^5
         -38934*c4^3*c5^4+83790*c4^3*c5^5)*c3^2
   -9*c4*c5^3*(-1+c4)*(-9471*c4*c5^3+6664*c4*c5^4-360*c4+1000*c4*c5+2532*c4*c5^2
         -180*c5-1764*c5^3+980*c5^4+1000*c5^2)*c3
   +9000*c4^3*c5^5-15876*c4^3*c5^6+8820*c4^3*c5^7-8820*c4^2*c5^7+1620*c4^2*c5^4
         -9000*c4^2*c5^5+15876*c4^2*c5^6-1620*c4^3*c5^4:

RK79_td1 := (c3,c4,c5,a52)-> 112*a52^2*(-10*c4+15*c4^2+2)*c3^6
   -8*a52*(24*a52+140*c4^2*a52-104*c4*a52+126*c4*c5-315*c4*c5^2-315*c4^2*c5+945*c4^2*c5^2)*c3^5
   +(48*a52^2+144*c5^2*a52-336*c5^3*a52-192*c4*a52^2+224*c4^2*a52^2+1134*c4^2*c5^2-5670*c4^2*c5^3
     +8505*c4^2*c5^4+5040*c4^2*c5^3*a52+2520*c4^2*c5^2*a52-1584*c4^2*c5*a52-840*c4*c5^3*a52
     +720*c4*c5*a52-1152*c4*c5^2*a52)*c3^4
   -6*c5*(315*c4*c5^4+1890*c4^2*c5^4-234*c4*c5^3-630*c4^2*c5^3-24*c5^2*a52-152*c4*c5^2*a52
     -216*c4^2*c5^2+54*c4*c5^2+560*c4^2*c5^2*a52+12*c5*a52+24*c4*c5*a52+108*c4^2*c5
     -96*c4^2*c5*a52-48*c4^2*a52+24*c4*a52)*c3^3
   +3*c5^2*(420*c4*c5^4+1260*c4^2*c5^4+42*c5^4+1260*c4^2*c5^3-144*c4*c5^3-36*c5^3-1224*c4^2*c5^2
     -72*c4*c5^2+9*c5^2-96*c4*c5*a52+224*c4^2*c5*a52+216*c4^2*c5+36*c4*c5+48*c4*a52
     -96*c4^2*a52+36*c4^2)*c3^2
   -36*c4*c5^3*(14*c5^3+70*c4*c5^3-12*c5^2-24*c4*c5^2+3*c5-12*c4*c5+6*c4)*c3
   +504*c4^2*c5^6+108*c4^2*c5^4-432*c4^2*c5^5:

RK79_L1 := (c3,c4,c5,a52)->RK79_tn1(c3,c4,c5,a52)/RK79_td1(c3,c4,c5,a52)/105:

RK79_L2 := (c3,c4,c5,a52)->(25200*c3^7*c4^3*c5*a52^2-73152*c3^4*c4^2*c5^2*a52-83700*c3^4*c4^2*c5^3*a52
   +139140*c3^4*c4^2*c5^4*a52-26256*c3^4*c4^3*c5*a52-2688*c3^4*c4^3*c5*a52^2+47268*c3^4*c4^3*c5^2*a52
   +85380*c3^4*c4^3*c5^3*a52-117600*c3^4*c4^3*c5^4*a52+19764*c3^5*c4*c5*a52+16432*c3^5*c4*c5*a52^2
   -74592*c3^5*c4*c5^2*a52+50580*c3^5*c4*c5^3*a52+11340*c3^5*c4*c5^4*a52-74100*c3^5*c4^2*c5*a52
   -31296*c3^5*c4^2*c5*a52^2+248688*c3^5*c4^2*c5^2*a52-111540*c3^5*c4^2*c5^3*a52-86100*c3^5*c4^2*c5^4*a52
   +59148*c3^5*c4^3*c5*a52+18480*c3^5*c4^3*c5*a52^2-188280*c3^5*c4^3*c5^2*a52+69300*c3^5*c4^3*c5^3*a52
   +75600*c3^5*c4^3*c5^4*a52-12096*c3^6*c4*c5*a52-31296*c3^6*c4*c5*a52^2+47628*c3^6*c4*c5^2*a52
   -39060*c3^6*c4*c5^3*a52+63360*c3^6*c4^2*c5*a52^2-181440*c3^6*c4^2*c5^2*a52+144900*c3^6*c4^2*c5^3*a52
   +4752*c3^2*c4*c5^3*a52+47628*c3^6*c4^2*c5*a52-39060*c3^6*c4^3*c5*a52-39200*c3^6*c4^3*c5*a52^2
   +144900*c3^6*c4^3*c5^2*a52-113400*c3^6*c4^3*c5^3*a52+18480*c3^7*c4*c5*a52^2-39200*c3^7*c4^2*c5*a52^2
   -1440*c3^2*c4*c5^2*a52-23544*c3^4*c4*c5^4*a52+34092*c3^4*c4^2*c5*a52+4784*c3^4*c4^2*c5*a52^2
   -1524*c3^4*c4*c5^3*a52-3528*c3^2*c4*c5^4*a52+4752*c3^2*c4^2*c5^2*a52-15360*c3^2*c4^2*c5^3*a52
   +11208*c3^2*c4^2*c5^4*a52-3528*c3^2*c4^3*c5^2*a52+11208*c3^2*c4^3*c5^3*a52-8064*c3^2*c4^3*c5^4*a52
   +1440*c3^3*c4*c5*a52+936*c3^3*c4*c5^2*a52-17472*c3^3*c4*c5^3*a52+16620*c3^3*c4*c5^4*a52
   -4752*c3^3*c4^2*c5*a52-11556*c3^3*c4^2*c5^2*a52+81924*c3^3*c4^2*c5^3*a52-70632*c3^3*c4^2*c5^4*a52
   +3528*c3^3*c4^3*c5*a52+11016*c3^3*c4^3*c5^2*a52-66600*c3^3*c4^3*c5^3*a52+55440*c3^3*c4^3*c5^4*a52
   -9720*c3^4*c4*c5*a52-2688*c3^4*c4*c5*a52^2+29916*c3^4*c4*c5^2*a52-270*c3^2*c5^4+1512*c3^2*c5^5
   -2691*c3^2*c5^6+1512*c3^2*c5^7-6048*c4^3*c5^5+10764*c4^3*c5^6-6048*c4^3*c5^7+2688*c3^5*a52^2
   -4784*c3^6*a52^2-10764*c4^2*c5^6+6048*c4^2*c5^7+1080*c4^3*c5^4-1080*c4^2*c5^4+6048*c4^2*c5^5
   +270*c3^3*c5^4-1512*c3^3*c5^5+2691*c3^3*c5^6-1512*c3^3*c5^7+2688*c3^7*a52^2-480*c3^4*a52^2
   -1080*c3^2*c4*c5^3-31851*c3^2*c4*c5^6+2754*c3^2*c4*c5^4-1080*c3^2*c4^2*c5^2+9234*c3^2*c4*c5^5
   +41580*c3*c4^3*c5^7-2160*c3*c4^3*c5^3-58320*c3*c4^3*c5^6+1080*c3^2*c4^3*c5^2+96822*c3^2*c4^3*c5^5
   -9720*c3^2*c4^2*c5^6+39870*c3^2*c4^3*c5^6+1764*c3^3*c5^4*a52-69084*c3^2*c4^3*c5^4+8694*c3^2*c4^3*c5^3
   +66150*c3^2*c4^2*c5^7+65232*c3^2*c4^2*c5^4-103734*c3^2*c4^2*c5^5+22302*c3^2*c4*c5^7
   -7452*c3^2*c4^2*c5^3-198450*c3^5*c4^3*c5^4-127575*c3^5*c4^2*c5^5-13608*c3^5*c4^3*c5^2
   +31296*c3^6*c4*a52^2+4784*c3^6*c5*a52^2+127575*c3^5*c4^3*c5^5-63360*c3^6*c4^2*a52^2+1080*c3*c4*c5^4
   +198450*c3^5*c4^2*c5^4+93555*c3^5*c4^3*c5^3+15444*c3*c4^3*c5^5-18480*c3^5*c4^3*a52^2-7128*c3*c4^2*c5^4
   +10764*c3*c4*c5^6+6048*c3*c4^3*c5^4-9396*c3*c4^2*c5^5+47556*c3*c4^2*c5^6-35532*c3*c4^2*c5^7
   +39200*c3^7*c4^2*a52^2-18480*c3^7*c4*a52^2+39200*c3^6*c4^3*a52^2-2688*c3^7*c5*a52^2-6048*c3*c4*c5^7
   +2160*c3*c4^2*c5^3-25200*c3^7*c4^3*a52^2-6048*c3*c4*c5^5+1764*c3^5*c5^2*a52-170100*c3^4*c4^3*c5^6
   +20682*c3^4*c4^3*c5^2+155250*c3^4*c4^3*c5^4-2688*c3^5*c5*a52^2+4032*c3^5*c5^4*a52+89775*c3^4*c4^3*c5^5
   -168966*c3^4*c4^2*c5^4+118827*c3^4*c4^2*c5^3-4784*c3^4*c4^2*a52^2-20250*c3^4*c4^2*c5^2
   +2688*c3^4*c4^3*a52^2-118665*c3^4*c4^3*c5^3-54675*c3^4*c4^2*c5^5+146475*c3^4*c4^2*c5^6
   +31296*c3^5*c4^2*a52^2+13608*c3^5*c4^2*c5^2-16432*c3^5*c4*a52^2-5604*c3^5*c5^3*a52-93555*c3^5*c4^2*c5^3
   -40950*c3^3*c4^2*c5^7+8532*c3^3*c4^2*c5^2-27513*c3^3*c4^2*c5^3-152235*c3^3*c4^2*c5^6
   +254529*c3^3*c4^2*c5^5-62010*c3^3*c4^2*c5^4-2376*c3^3*c5^3*a52+4806*c3^3*c4*c5^3-88200*c3^2*c4^3*c5^7
   +720*c3^3*c5^2*a52-1620*c3^3*c4*c5^6-18900*c3^3*c4*c5^7-26784*c3^3*c4*c5^4+39744*c3^3*c4*c5^5
   +2688*c3^4*c4*a52^2+7680*c3^4*c5^3*a52-2376*c3^4*c5^2*a52-5604*c3^4*c5^4*a52-48681*c3^4*c4*c5^5
   +25002*c3^4*c4*c5^4-3969*c3^4*c4*c5^3+29295*c3^4*c4*c5^6+23706*c3^3*c4^3*c5^3+87912*c3^3*c4^3*c5^4
   -8694*c3^3*c4^3*c5^2+160650*c3^3*c4^3*c5^6+480*c3^4*c5*a52^2-297810*c3^3*c4^3*c5^5+56700*c3^3*c4^3*c5^7)
   /(2520*c3^4*c4^2*c5^2*a52+5040*c3^4*c4^2*c5^3*a52-1008*c3^5*c4*c5*a52+2520*c3^5*c4*c5^2*a52
   +2520*c3^5*c4^2*c5*a52-7560*c3^5*c4^2*c5^2*a52-288*c3^2*c4*c5^3*a52+144*c3^2*c4*c5^2*a52
   -1584*c3^4*c4^2*c5*a52-840*c3^4*c4*c5^3*a52-288*c3^2*c4^2*c5^2*a52+672*c3^2*c4^2*c5^3*a52
   -144*c3^3*c4*c5*a52-144*c3^3*c4*c5^2*a52+912*c3^3*c4*c5^3*a52+288*c3^3*c4^2*c5*a52
   +576*c3^3*c4^2*c5^2*a52-3360*c3^3*c4^2*c5^3*a52+720*c3^4*c4*c5*a52-1152*c3^4*c4*c5^2*a52
   +27*c3^2*c5^4-108*c3^2*c5^5+126*c3^2*c5^6-192*c3^5*a52^2+224*c3^6*a52^2+504*c4^2*c5^6+108*c4^2*c5^4
   -432*c4^2*c5^5+48*c3^4*a52^2+108*c3^2*c4*c5^3+1260*c3^2*c4*c5^6-216*c3^2*c4*c5^4+108*c3^2*c4^2*c5^2
   -432*c3^2*c4*c5^5+3780*c3^2*c4^2*c5^6-3672*c3^2*c4^2*c5^4+3780*c3^2*c4^2*c5^5+648*c3^2*c4^2*c5^3
   -1120*c3^6*c4*a52^2+1680*c3^6*c4^2*a52^2-108*c3*c4*c5^4+432*c3*c4^2*c5^4-504*c3*c4*c5^6
   +864*c3*c4^2*c5^5-2520*c3*c4^2*c5^6-216*c3*c4^2*c5^3+432*c3*c4*c5^5+8505*c3^4*c4^2*c5^4
   -5670*c3^4*c4^2*c5^3+224*c3^4*c4^2*a52^2+1134*c3^4*c4^2*c5^2-1120*c3^5*c4^2*a52^2+832*c3^5*c4*a52^2
   -648*c3^3*c4^2*c5^2+1296*c3^3*c4^2*c5^3-11340*c3^3*c4^2*c5^5+3780*c3^3*c4^2*c5^4+144*c3^3*c5^3*a52
   -324*c3^3*c4*c5^3-72*c3^3*c5^2*a52+1404*c3^3*c4*c5^4-1890*c3^3*c4*c5^5-192*c3^4*c4*a52^2
   -336*c3^4*c5^3*a52+144*c3^4*c5^2*a52)/60:

RK79_L4 := (c3,c4,c5,a52)->
   -1/60*(-7344*c4^3*c5^5+13248*c4^3*c5^6-7560*c4^3*c5^7+7560*c4^2*c5^7+1296*c4^3*c5^4-1296*c4^2*c5^4
   +324*c3^3*c5^4-1836*c3^3*c5^5+3312*c3^3*c5^6-1890*c3^3*c5^7-576*c3^4*a52^2+3360*c3^7*a52^2
   -5888*c3^6*a52^2+7344*c4^2*c5^5-13248*c4^2*c5^6+3264*c3^5*a52^2+1836*c3^2*c5^5-3312*c3^2*c5^6
   +1890*c3^2*c5^7-324*c3^2*c5^4-24570*c3^3*c4*c5^7-32940*c3^3*c4*c5^4+49626*c3^3*c4*c5^5
   +109890*c3^3*c4^3*c5^4-78075*c3^3*c4^2*c5^4+10368*c3^3*c4^2*c5^2-33804*c3^3*c4^2*c5^3
   -114660*c3^2*c4^3*c5^7+864*c3^3*c5^2*a52+170100*c3^5*c4^3*c5^5-85230*c3^2*c4^3*c5^4
   -2025*c3^3*c4*c5^6+5832*c3^3*c4*c5^3-2880*c3^3*c5^3*a52+2160*c3^3*c5^4*a52+3264*c3^4*c4*a52^2
   +576*c3^4*c5*a52^2+209790*c3^3*c4^3*c5^6+9408*c3^4*c5^3*a52+31077*c3^4*c4*c5^4-61560*c3^4*c4*c5^5
   -2880*c3^4*c5^2*a52-6936*c3^4*c5^4*a52-10584*c3^3*c4^3*c5^2-197190*c3^3*c4^2*c5^6
   -56070*c3^3*c4^2*c5^7+324810*c3^3*c4^2*c5^5+75600*c3^3*c4^3*c5^7-379890*c3^3*c4^3*c5^5
   +29322*c3^3*c4^3*c5^3+50960*c3^6*c4^3*a52^2-3360*c3^7*c5*a52^2+39120*c3^6*c4*a52^2
   -80640*c3^6*c4^2*a52^2-33600*c3^7*c4^3*a52^2+52920*c3*c4^3*c5^7-23520*c3^7*c4*a52^2
   +50960*c3^7*c4^2*a52^2-170100*c3^5*c4^2*c5^5-23520*c3^5*c4^3*a52^2+119070*c3^5*c4^3*c5^3
   +257985*c3^5*c4^2*c5^4-17010*c3^5*c4^3*c5^2-257985*c3^5*c4^3*c5^4+5888*c3^6*c5*a52^2
   -8640*c3*c4^2*c5^4-2592*c3*c4^3*c5^3-45360*c3*c4^2*c5^7+19008*c3*c4^3*c5^5+39120*c3^5*c4^2*a52^2
   -119070*c3^5*c4^2*c5^3+7344*c3*c4^3*c5^4+17010*c3^5*c4^2*c5^2-7560*c3*c4*c5^7+59652*c3*c4^2*c5^6
   +2592*c3*c4^2*c5^3+1296*c3*c4*c5^4+13248*c3*c4*c5^6+10584*c3^2*c4^3*c5^3-13500*c3^2*c4^2*c5^6
   +50400*c3^2*c4^3*c5^6-129384*c3^2*c4^2*c5^5+86940*c3^2*c4^2*c5^7-7344*c3*c4*c5^5
   +1296*c3^2*c4^3*c5^2+121500*c3^2*c4^3*c5^5-11664*c3*c4^2*c5^5+25542*c3^4*c4^3*c5^2
   -149040*c3^4*c4^3*c5^3+197505*c3^4*c4^2*c5^6+199395*c3^4*c4^3*c5^4-226800*c3^4*c4^3*c5^6
   -6936*c3^5*c5^3*a52+115290*c3^4*c4^3*c5^5-5888*c3^4*c4^2*a52^2+148554*c3^4*c4^2*c5^3
   -4860*c3^4*c4*c5^3+37800*c3^4*c4*c5^6-214245*c3^4*c4^2*c5^4-73710*c3^4*c4^2*c5^5
   -24948*c3^4*c4^2*c5^2+3360*c3^4*c4^3*a52^2-39762*c3^2*c4*c5^6+11340*c3^2*c4*c5^5-1296*c3^2*c4*c5^3
   +3348*c3^2*c4*c5^4+80352*c3^2*c4^2*c5^4-9072*c3^2*c4^2*c5^3+28350*c3^2*c4*c5^7-1296*c3^2*c4^2*c5^2
   +5040*c3^5*c5^4*a52-20224*c3^5*c4*a52^2+2160*c3^5*c5^2*a52-3264*c3^5*c5*a52^2-72900*c3*c4^3*c5^6
   +33600*c3^7*c4^3*c5*a52^2-10080*c3^2*c4^3*c5^4*a52+1728*c3^3*c4*c5*a52+1152*c3^3*c4*c5^2*a52
   +5760*c3^2*c4^2*c5^2*a52-18816*c3^2*c4^2*c5^3*a52+13872*c3^2*c4^2*c5^4*a52-4320*c3^2*c4^3*c5^2*a52
   +13872*c3^2*c4^3*c5^3*a52-21504*c3^3*c4*c5^3*a52+20688*c3^3*c4*c5^4*a52-15120*c3^6*c4*c5*a52
   -39120*c3^6*c4*c5*a52^2+60480*c3^6*c4*c5^2*a52-50400*c3^6*c4*c5^3*a52+60480*c3^6*c4^2*c5*a52
   +80640*c3^6*c4^2*c5*a52^2-234360*c3^6*c4^2*c5^2*a52+190260*c3^6*c4^2*c5^3*a52
   -50400*c3^6*c4^3*c5*a52-50960*c3^6*c4^3*c5*a52^2+190260*c3^6*c4^3*c5^2*a52
   -151200*c3^6*c4^3*c5^3*a52+23520*c3^7*c4*c5*a52^2-50960*c3^7*c4^2*c5*a52^2-1728*c3^2*c4*c5^2*a52
   +5760*c3^2*c4*c5^3*a52-4320*c3^2*c4*c5^4*a52+100800*c3^5*c4^3*c5^4*a52-241920*c3^5*c4^3*c5^2*a52
   +89460*c3^5*c4^3*c5^3*a52-1560*c3^4*c4*c5^3*a52-30060*c3^4*c4*c5^4*a52+41904*c3^4*c4^2*c5*a52
   +5888*c3^4*c4^2*c5*a52^2-90864*c3^4*c4^2*c5^2*a52-106740*c3^4*c4^2*c5^3*a52
   +178920*c3^4*c4^2*c5^4*a52-32664*c3^4*c4^3*c5*a52-3360*c3^4*c4^3*c5*a52^2+59400*c3^4*c4^3*c5^2*a52
   +110040*c3^4*c4^3*c5^3*a52-152880*c3^4*c4^3*c5^4*a52+24336*c3^5*c4*c5*a52+20224*c3^5*c4*c5*a52^2
   -93024*c3^5*c4*c5^2*a52+63540*c3^5*c4*c5^3*a52+15120*c3^5*c4*c5^4*a52-92496*c3^5*c4^2*c5*a52
   -39120*c3^5*c4^2*c5*a52^2+314820*c3^5*c4^2*c5^2*a52-141960*c3^5*c4^2*c5^3*a52
   -113820*c3^5*c4^2*c5^4*a52+74880*c3^5*c4^3*c5*a52+23520*c3^5*c4^3*c5*a52^2-5760*c3^3*c4^2*c5*a52
   -14256*c3^3*c4^2*c5^2*a52+101976*c3^3*c4^2*c5^3*a52-88920*c3^3*c4^2*c5^4*a52+4320*c3^3*c4^3*c5*a52
   +13752*c3^3*c4^3*c5^2*a52-83880*c3^3*c4^3*c5^3*a52+70560*c3^3*c4^3*c5^4*a52-11808*c3^4*c4*c5*a52
   -3264*c3^4*c4*c5*a52^2+36720*c3^4*c4*c5^2*a52)/(108*c4^2*c5^4+48*c3^4*a52^2+224*c3^6*a52^2
   -432*c4^2*c5^5+504*c4^2*c5^6-192*c3^5*a52^2-108*c3^2*c5^5+126*c3^2*c5^6+27*c3^2*c5^4
   +1404*c3^3*c4*c5^4-1890*c3^3*c4*c5^5+3780*c3^3*c4^2*c5^4-648*c3^3*c4^2*c5^2+1296*c3^3*c4^2*c5^3
   -72*c3^3*c5^2*a52-324*c3^3*c4*c5^3+144*c3^3*c5^3*a52-192*c3^4*c4*a52^2-336*c3^4*c5^3*a52
   +144*c3^4*c5^2*a52-11340*c3^3*c4^2*c5^5-1120*c3^6*c4*a52^2+1680*c3^6*c4^2*a52^2+432*c3*c4^2*c5^4
   -1120*c3^5*c4^2*a52^2-2520*c3*c4^2*c5^6-216*c3*c4^2*c5^3-108*c3*c4*c5^4-504*c3*c4*c5^6
   +3780*c3^2*c4^2*c5^6+3780*c3^2*c4^2*c5^5+432*c3*c4*c5^5+864*c3*c4^2*c5^5+224*c3^4*c4^2*a52^2
   -5670*c3^4*c4^2*c5^3+8505*c3^4*c4^2*c5^4+1134*c3^4*c4^2*c5^2+1260*c3^2*c4*c5^6-432*c3^2*c4*c5^5
   +108*c3^2*c4*c5^3-216*c3^2*c4*c5^4-3672*c3^2*c4^2*c5^4+648*c3^2*c4^2*c5^3+108*c3^2*c4^2*c5^2
   +832*c3^5*c4*a52^2-144*c3^3*c4*c5*a52-144*c3^3*c4*c5^2*a52-288*c3^2*c4^2*c5^2*a52+
   672*c3^2*c4^2*c5^3*a52+912*c3^3*c4*c5^3*a52+144*c3^2*c4*c5^2*a52-288*c3^2*c4*c5^3*a52
   -840*c3^4*c4*c5^3*a52-1584*c3^4*c4^2*c5*a52+2520*c3^4*c4^2*c5^2*a52+5040*c3^4*c4^2*c5^3*a52
   -1008*c3^5*c4*c5*a52+2520*c3^5*c4*c5^2*a52+2520*c3^5*c4^2*c5*a52-7560*c3^5*c4^2*c5^2*a52
   +288*c3^3*c4^2*c5*a52+576*c3^3*c4^2*c5^2*a52-3360*c3^3*c4^2*c5^3*a52+720*c3^4*c4*c5*a52
   -1152*c3^4*c4*c5^2*a52):

RK79_b9 := (c3,c4,c5,a52,c7,c8) -> 
      (RK79_L1(c3,c4,c5,a52) + RK79_L2(c3,c4,c5,a52)*(c7+c8) + RK79_L4(c3,c4,c5,a52)*c7*c8)
     /((c8-1)*(c7-1)*(RK79_c6(c3,c4,c5,a52)-1)*(c5-1)*(c4-1)*(c3-1));

# ---------------------- car1 ---------------------------------------------------------------
# RK79_car1 (c3,c4,c5,a52,c7,c8) -> [b3,b4,b5,b6,b7,b8,b9]
RK79_car1 := proc(c3,c4,c5,c6,c7,c8,b9) local b3,b4,b5,b6,b7,b8, solb, 
     car1,car2,car3,car4,car5,car6;
  car1:= b3*c3   + b4*c4   + b5*c5   + b6*c6   + b7*c7   + b8*c8   + b9 - 1/2;
  car2:= b3*c3^2 + b4*c4^2 + b5*c5^2 + b6*c6^2 + b7*c7^2 + b8*c8^2 + b9 - 1/3;
  car3:= b3*c3^3 + b4*c4^3 + b5*c5^3 + b6*c6^3 + b7*c7^3 + b8*c8^3 + b9 - 1/4;
  car4:= b3*c3^4 + b4*c4^4 + b5*c5^4 + b6*c6^4 + b7*c7^4 + b8*c8^4 + b9 - 1/5;
  car5:= b3*c3^5 + b4*c4^5 + b5*c5^5 + b6*c6^5 + b7*c7^5 + b8*c8^5 + b9 - 1/6;
  car6:= b3*c3^6 + b4*c4^6 + b5*c5^6 + b6*c6^6 + b7*c7^6 + b8*c8^6 + b9 - 1/7;
#   ------------ b[i] -----------------------------
  solb:=solve({car1, car2, car3, car4,  car5,  car6},
              {b3,   b4,   b5,   b6,    b7,    b8});
  assign(solb);
  [b3,b4,b5,b6,b7,b8];
end proc:

# ---------------------- car1 ---------------------------------------------------------------
RK79_c7c8 := proc(c3,c4,c5,a52) local c6,c7,c8,b3,b4,b5,b6,b7,b8,b9, solb,
             a32,a42,a62,a72,a82,a92, eq0, eq1, eq2, eq3, eq17, sole; 
  c6 := RK79_c6  (c3,c4,c5,a52);
  #print(`  c6= `, c6);
  b9 := RK79_b9  (c3,c4,c5,a52,c7,c8); 
  #print(`  b9= `, b9);
  solb := RK79_car1(c3,c4,c5,c6,c7,c8,b9);
  b3 := solb[1]:
  b4 := solb[2]:
  b5 := solb[3]:
  b6 := solb[4]:
  b7 := solb[5]:
  b8 := solb[6]:
  #print(`  b8= `, b8);
  a32:= 3/4*c3;
  a42:= 3*c4^2*(3*c3-2*c4)/c3^2/4;
  #c7 := 7:    #  !!!!!!!!11
  #c8 := 8:    #  !!!!!!!!11
# eq1,..eq5
  eq0 := b3*a32      + b4*a42      + b5*a52      + b6*a62      + b7*a72      + b8*a82      + b9*a92;
  eq1 := b3*a32*c3   + b4*a42*c4   + b5*a52*c5   + b6*a62*c6   + b7*a72*c7   + b8*a82*c8   + b9*a92;
  eq2 := b3*a32*c3^2 + b4*a42*c4^2 + b5*a52*c5^2 + b6*a62*c6^2 + b7*a72*c7^2 + b8*a82*c8^2 + b9*a92;
  eq3 := b3*a32*c3^3 + b4*a42*c4^3 + b5*a52*c5^3 + b6*a62*c6^3 + b7*a72*c7^3 + b8*a82*c8^3 + b9*a92;
  eq17:= b3*a32^2    + b4*a42^2    + b5*a52^2    + b6*a62^2    + b7*a72^2    + b8*a82^2    + b9*a92^2;

  eq0 := factor(eval(eq0));
  eq1 := factor(eval(eq1));
  eq2 := factor(eval(eq2));
  eq3 := factor(eval(eq3));
  eq17 := factor(eval(eq17));

  #print( 770, eq0 );
  #print( 771, eq1 );
  #print( 772, eq2 );
  #print( 773, eq3 );
  #print( 7717,eq17);
 
  sole:=solve({eq0, eq1, eq2, eq3},
              {a62, a72, a82, a92}); #print(sole);
  assign(sole);
  #print(`  a92= `, a92);
  eq17:=factor(eval(eq17));
  eq17:=numer(eq17);
  #print(`  eq17= `, eq17);
  unapply(eq17,c7,c8);
end proc:

# Rational point on biquadric curve

#  point P(x,y) -> point(x1,y)
nextX:=proc(f,P) local x,y, fx0, x1;
  if simplify( f(P[1],P[2]) ) <> 0 then ERROR(`nextX: f(x,y) <>0`); end if;
  fx0:=factor(f(x,P[2])/(x-P[1])); #print(fx0);
  x1 := solve(fx0,x); #print(x1);
  [x1, P[2]];
end proc:

#  point P(x,y) -> point(x,y1)
nextY:=proc(f,P) local x,y, fy0, y1;
  if simplify( f(P[1],P[2]) ) <> 0 then ERROR(`nextX: f(x,y) <>0`); end if;
  fy0:=factor(f(P[1],y)/(y-P[2])); #print(fy0);
  y1 := solve(fy0,y); #print(y1);
  [P[1], y1];
end proc:

# points P1(x1,y1),P2(x2,y2) -> point P3(x3,y3)
point3:=proc(f,P1,P2) local g, x1,y1, x2,y2,  gt,t, x,y;
  #print(777, evalf(P1), evalf(P2));
  x1:=P1[1]; y1:=P1[2]; 
  x2:=P2[1]; y2:=P2[2]; 
  #print(778, x1, y1, x2, y2);

  if simplify( f(x1,y1)) <> 0 then ERROR(`point3: f(x1,y1) <>0`); end if;
  if simplify( f(x2,y2)) <> 0 then ERROR(`point3: f(x2,y2) <>0`); end if;
  gt:= f(x1 + t*(x2-x1), y1 + t*(y2-y1));
  gt:= factor(gt/t/(t-1));
  #print(`  gt= (linear of t)=  `, evalf(gt));
  #print(779, degree(gt,t));
  if degree(gt,t) <> 1 then return [0,0]; end if;
  t:=solve(gt,t);
  #print(`  t= `, t, evalf(t));
  #print(x1 + t*(x2-x1), y1 + t*(y2-y1), evalf(x1 + t*(x2-x1)), evalf(y1 + t*(y2-y1)));
  [x1 + t*(x2-x1), y1 + t*(y2-y1)];
end proc;

# about 4N rational point on the curve f(x,y)=0
ratPoints:=proc(f,N1::integer) local N, P, Q, x, y, sol, x0, y0, k, f3, f3xy,  P2, qt, qx, qy, i,j, p0;
  N:=N1;
  if type(N,odd) then N:=N+1; end if;
  P:=Array(0..N);
  Q:=Array(1..N);
  sol:=solve(f(1,y));
  if nops([sol]) <> 2 then ERROR(`ratPoint: f(1,y) not solved`); end if; 
  x0 := sol[1]; y0 := 1;
  P[0] := [x0, y0];
  for k from 1 to N by 2 do
    P[k  ] := nextX(f, P[k-1]): #print(k  , evalf(P[k  ]), length(P[k  ]) );
    P[k+1] := nextY(f, P[k  ]): #print(k+1, evalf(P[k+1]), length(P[k+1]) );
  end do:
  #for k from 0 to N do print(k, f(P[k][1], P[k][2])); end do;

  f3xy := factor( f(x+x0, x/y+y0)*y^2/x); # cubic polynom(x,y)
  f3 := unapply(f3xy,x,y);
  #print(703, f3xy);
  #print(704, f3(x,y));

  for k from 2 to N do
  #           u = x-x0,    v = x/(y-y0)
    Q[k] := [ P[k][1]-x0, (P[k][1]-x0)/(P[k][2]-y0)];
  end do:
  #for k from 2 to N do print(703, k , f3(Q[k][1], Q[k][2])); end do;

  P2:=[]:
  for k from 2 to N do
    if P[k][1]=0 or P[k][1]=1 or P[k][2]=0 or P[k][2]=1 then next end if;
    P2 := [op(P2), P[k]];
  end do:

  for i from 2 to N do 
    for j from 2 to i-1 do 
      qt:=point3(f3, Q[i], Q[j]);
      qx:=qt[1]: qy:=qt[2]:
      if qx=0 or qx=1 or qy=0 or qy=1 then next; end if;
      p0:=[qx+x0, qx/qy+y0];
      if p0 in P2 then next; end if;
      P2 := [op(P2), p0];
    end do:
  end do:

  P2;
end proc:
# ==========================================================
RK_79a := proc(c3,c4,c5,a52,c7,c8)::Matrix; local s, a,b,c,A, eq0,eq1,eq2,eq3,eq17,sole,
      solb, i,j, e26, e27, e37_3, e37_4, e37_5, e37_0, sol, 
      eq, w4s7,w4s8,w4s9, w5, k, Dnx, x4, x2, x4a, x3, Zt;

  s:=9; 
  c[3]:=c3: c[4]:=c4: c[5]:=c5: c[7]:=c7: c[8]:=c8:
  a[5,2]:=a52:

  c[6] := RK79_c6  (c3,c4,c5,a52);
  b[9] := RK79_b9  (c3,c4,c5,a52,c7,c8); 
  solb := RK79_car1(c3,c4,c5,c[6],c7,c8,b[9]);
  b[3] := solb[1]:
  b[4] := solb[2]:
  b[5] := solb[3]:
  b[6] := solb[4]:
  b[7] := solb[5]:
  b[8] := solb[6]:

  c[1]:=0: c[s]:=1: c[s+1]:=1: b[2]:=0:
  c[2]  :=2/3*c3: 
  a[3,2]:= 3/4*c3;
  a[4,2]:= 3*c4^2*(3*c3-2*c4)/c3^2/4;
  a[4,3] := c4^2*(c4-c3)/c3^2;
  # eq1,..eq5
  eq0 := b[3]*a[3,2]        + b[4]*a[4,2]        + b[5]*a[5,2]        + b[6]*a[6,2]        + b[7]*a[7,2]        + b[8]*a[8,2]        + b[9]*a[9,2];
  eq1 := b[3]*a[3,2]*c[3]   + b[4]*a[4,2]*c[4]   + b[5]*a[5,2]*c[5]   + b[6]*a[6,2]*c[6]   + b[7]*a[7,2]*c[7]   + b[8]*a[8,2]*c[8]   + b[9]*a[9,2];
  eq2 := b[3]*a[3,2]*c[3]^2 + b[4]*a[4,2]*c[4]^2 + b[5]*a[5,2]*c[5]^2 + b[6]*a[6,2]*c[6]^2 + b[7]*a[7,2]*c[7]^2 + b[8]*a[8,2]*c[8]^2 + b[9]*a[9,2];
  eq3 := b[3]*a[3,2]*c[3]^3 + b[4]*a[4,2]*c[4]^3 + b[5]*a[5,2]*c[5]^3 + b[6]*a[6,2]*c[6]^3 + b[7]*a[7,2]*c[7]^3 + b[8]*a[8,2]*c[8]^3 + b[9]*a[9,2];
  eq17:= b[3]*a[3,2]^2      + b[4]*a[4,2]^2      + b[5]*a[5,2]^2      + b[6]*a[6,2]^2      + b[7]*a[7,2]^2      + b[8]*a[8,2]^2      + b[9]*a[9,2]^2;

  eq0 := factor(eval(eq0 ));
  eq1 := factor(eval(eq1 ));
  eq2 := factor(eval(eq2 ));
  eq3 := factor(eval(eq3 ));
  eq17:= factor(eval(eq17));

  #print( 770, eq0 );
  #print( 771, eq1 );
  #print( 772, eq2 );
  #print( 773, eq3 );
  #print( 7717,eq17);
 
  sole:=solve({eq0,    eq1,    eq2,    eq3},
              {a[6,2], a[7,2], a[8,2], a[9,2]}); #print(sole);
  assign(sole);
  #if factor(eval(eq17)) <> 0 then ERROR(`RK_79a: equation RK79_c7c8 not holds`); end if;

  # ---------- b[i], c[i], a[i,2] are ready -----------------------

  a[5,3] := (9*c[5]^2*c[4]-12*a[5,2]*c[3]*c[4]+8*a[5,2]*c[3]^2-6*c[5]^3)/(c[4]-c[3])/c[3]/18;
  a[5,4] := (6*c[5]^3-9*c[5]^2*c[3]+4*a[5,2]*c[3]^2)/c[4]/(c[4]-c[3])/18;

  # last equations in 1st simplifying assumption:
  a[9,8] := (1-c[8])*b[8]/b[9];

  #-------------------------------------------------------------------------------------------------
  e26 := 1/2*c[6]^2-2/3*a[6,2]*c[3]-a[6,3]*c[3]-a[6,4]*c[4]-a[6,5]*c[5];
  e27 := c[6]^3-c[6]^2*c[3]-a[6,3]*c[3]^2-3*a[6,4]*c[4]^2+2*a[6,4]*c[4]*c[3]-3*a[6,5]*c[5]^2+2*a[6,5]*c[5]*c[3];
  e37_3:= -9*c[3]^3*(12*a[5,2]*c[3]^3*c[4]+12*a[5,2]*c[3]^3*c[5]+4*a[5,2]*c[3]^2*c[5]*c[4]
          -16*c[4]^2*a[5,2]*c[3]^2-54*c[3]^2*c[5]^2*c[4]+36*c[4]*c[3]*c[5]^3+36*c[4]^2*c[5]^2*c[3]
          -9*c[5]^4*c[3]-24*c[4]^2*c[5]^3+6*c[5]^4*c[4]);
  
  e37_4:= -9*c[4]^2*(-24*c[4]^2*a[5,2]*c[3]^3-12*c[5]^4*c[4]^2+36*c[4]^2*c[5]^2*c[3]^2
          -8*a[5,2]*c[3]^2*c[5]*c[4]^2-54*c[3]^3*c[5]^2*c[4]+20*c[3]^4*a[5,2]*c[4]+36*a[5,2]*c[3]^3*c[5]*c[4]
          +36*c[3]*c[5]^4*c[4]-24*c[4]*c[3]^2*c[5]^3+36*c[3]^3*c[5]^3-27*c[5]^4*c[3]^2-12*a[5,2]*c[3]^4*c[5]);
  
  e37_5:= -36*a[5,2]*c[3]^3*c[5]^4+324*c[3]^2*c[5]^5*c[4]-288*c[4]*c[3]*c[5]^6-216*c[4]^2*c[5]^5*c[3]
          +64*a[5,2]^2*c[3]^4*c[4]^2-80*a[5,2]^2*c[3]^5*c[4]+48*a[5,2]^2*c[3]^5*c[5]-144*a[5,2]*c[3]^4*c[5]^3
          +108*c[3]^2*c[5]^6+144*c[4]^2*c[5]^6-96*a[5,2]*c[3]^3*c[4]*c[5]^3+192*c[4]^2*a[5,2]*c[3]^2*c[5]^3
          -48*a[5,2]^2*c[3]^4*c[5]*c[4]+216*a[5,2]*c[3]^4*c[5]^2*c[4]-144*a[5,2]*c[3]^3*c[4]^2*c[5]^2
          -72*a[5,2]*c[3]^2*c[5]^4*c[4];
  
  e37_0:=-48*a[6,2]*c[3]^6*a[5,2]*c[4]-48*a[6,2]*c[3]^6*a[5,2]*c[5]+48*a[6,2]*c[3]^5*c[4]^2*a[5,2]
         +216*a[6,2]*c[3]^5*c[5]^2*c[4]-180*a[6,2]*c[3]^4*c[4]*c[5]^3-108*a[6,2]*c[3]^4*c[4]^2*c[5]^2
         +72*a[6,2]*c[3]^3*c[4]^2*c[5]^3+12*a[5,2]*c[3]^3*c[4]*c[6]^4+36*a[5,2]*c[3]^3*c[5]*c[6]^4
         -24*c[4]^2*a[5,2]*c[3]^2*c[6]^4-81*c[3]^2*c[5]^2*c[4]*c[6]^4+72*c[4]*c[3]*c[5]^3*c[6]^4
         +54*c[4]^2*c[5]^2*c[3]*c[6]^4+36*a[6,2]*c[6]*c[3]^4*c[5]^3-36*a[6,2]*c[6]*c[3]^3*c[5]^4
         +36*a[6,2]*c[3]^4*c[5]^4-27*c[3]^2*c[5]^3*c[6]^4-36*c[4]^2*c[5]^3*c[6]^4
         -16*a[6,2]*c[6]*c[3]^4*c[4]^2*a[5,2]+32*a[6,2]*c[6]*c[3]^5*a[5,2]*c[4]
         -108*a[6,2]*c[6]*c[3]^4*c[5]^2*c[4]+36*a[6,2]*c[6]*c[3]^3*c[4]^2*c[5]^2
         +84*a[6,2]*c[6]*c[3]^3*c[4]*c[5]^3-24*a[6,2]*c[6]*c[3]^2*c[4]^2*c[5]^3;
  
  sol:=solve({e26,e27, e37_0+e37_3*a[6,3]+e37_4*a[6,4]+e37_5*a[6,5]},{a[6,3],a[6,4],a[6,5]});
  assign(sol);

  #-------------------------------------------------------------------------------------------------
  #e18:=add(w4s[i]*b[i]       , i=3..s);
  #e19:=add(w4s[i]*b[i]*c[i]  , i=3..s);
  #e20:=add(w4s[i]*b[i]*c[i]^2, i=3..s);

  eq[18]:= -b[3]*c[3]^4/9 - c[4]^2*(12*c[3]^2-20*c[4]*c[3]+9*c[4]^2)*b[4]/9
        +b[5]*(32/27*a[5,2]*c[3]^3+4*a[5,3]*c[3]^3+4*a[5,4]*c[4]^3-c[5]^4)
        +b[6]*(32/27*a[6,2]*c[3]^3+4*a[6,3]*c[3]^3+4*a[6,4]*c[4]^3+4*a[6,5]*c[5]^3-c[6]^4)
        +w4s7*b[7]+w4s8*b[8]+w4s9*b[9];
  eq[19]:= -b[3]*c[3]^5/9-c[4]^3*(12*c[3]^2-20*c[4]*c[3]+9*c[4]^2)*b[4]/9
        +(32/27*a[5,2]*c[3]^3+4*a[5,3]*c[3]^3+4*a[5,4]*c[4]^3-c[5]^4)*b[5]*c[5]
        +(32/27*a[6,2]*c[3]^3+4*a[6,3]*c[3]^3+4*a[6,4]*c[4]^3+4*a[6,5]*c[5]^3-c[6]^4)*b[6]*c[6]
        +w4s7*b[7]*c[7]+w4s8*b[8]*c[8]+w4s9*b[9];
  eq[20]:= -b[3]*c[3]^6/9-c[4]^4*(12*c[3]^2-20*c[4]*c[3]+9*c[4]^2)*b[4]/9
        +(32/27*a[5,2]*c[3]^3+4*a[5,3]*c[3]^3+4*a[5,4]*c[4]^3-c[5]^4)*b[5]*c[5]^2
        +(32/27*a[6,2]*c[3]^3+4*a[6,3]*c[3]^3+4*a[6,4]*c[4]^3+4*a[6,5]*c[5]^3-c[6]^4)*b[6]*c[6]^2
        +w4s7*b[7]*c[7]^2+w4s8*b[8]*c[8]^2+w4s9*b[9];
  sol:=solve({eq[18],eq[19],eq[20]},{w4s7,w4s8,w4s9});
  assign(sol);

  # ---------------------------------------------
  a[9,7] :=((1-c[7])*b[7] -b[8]*a[8,7])/b[9];   # C1 + C2*a[8,7]

  #-------------------------------------------------------------------------------------------------
  #     -w4s7 + w4[7]
  #     -w4s8 + w4[8]
  #     -w4s9 + w4[9]
  eq[21]:=-w4s7 + 4*(a[7,2]*c[2]^3+a[7,3]*c[3]^3+a[7,4]*c[4]^3+a[7,5]*c[5]^3+a[7,6]*c[6]^3)-c[7]^4;
  eq[22]:=-w4s8 + 4*(a[8,2]*c[2]^3+a[8,3]*c[3]^3+a[8,4]*c[4]^3+a[8,5]*c[5]^3+a[8,6]*c[6]^3+a[8,7]*c[7]^3)-c[8]^4;
  eq[23]:=-w4s9 + 4*(a[9,2]*c[2]^3+a[9,3]*c[3]^3+a[9,4]*c[4]^3+a[9,5]*c[5]^3+a[9,6]*c[6]^3+a[9,7]*c[7]^3+a[9,8]*c[8]^3)-1;

  # a[7,3], a[7,4]:
  eq[28]:= 1/2*c[7]^2-a[7,2]*c[2]-a[7,3]*c[3]-a[7,4]*c[4]-a[7,5]*c[5]-a[7,6]*c[6];
  eq[29]:= c[7]^2*(c[7]-c[3])-a[7,3]*c[3]^2-a[7,4]*c[4]*(3*c[4]-2*c[3])
                             -a[7,5]*c[5]*(3*c[5]-2*c[3])-a[7,6]*c[6]*(3*c[6]-2*c[3]);

  # a[8,3], a[8,4]:
  eq[30]:= 1/2*c[8]^2-a[8,2]*c[2]-a[8,3]*c[3]-a[8,4]*c[4]-a[8,5]*c[5]-a[8,6]*c[6]-a[8,7]*c[7];
  eq[31]:= c[8]^2*(c[8]-c[3])-a[8,3]*c[3]^2-a[8,4]*c[4]*(3*c[4]-2*c[3])-a[8,5]*c[5]*(3*c[5]-2*c[3])
                             -a[8,6]*c[6]*(3*c[6]-2*c[3])-a[8,7]*c[7]*(3*c[7]-2*c[3]);

  # a[9,3], a[9,4]:
  eq[32]:= 1/2-a[9,2]*c[2]-a[9,3]*c[3]-a[9,4]*c[4]-a[9,5]*c[5]-a[9,6]*c[6]-a[9,7]*c[7]-a[9,8]*c[8];
  eq[33]:= 1-c[3]-a[9,3]*c[3]^2-a[9,4]*c[4]*(3*c[4]-2*c[3])-a[9,5]*c[5]*(3*c[5]-2*c[3])-a[9,6]*c[6]*(3*c[6]-2*c[3])
                 -a[9,7]*c[7]*(3*c[7]-2*c[3])-a[9,8]*c[8]*(3*c[8]-2*c[3]);

  #s1:=s+1;
  w5:= Array(1..s+1);
  w5[2] := -c[2]^5;
  for k from 3 to s+1 do 
    w5[k] := 5*add( a[k,i]*c[i]^4,i=2..k-1 ) - c[k]^5;
  end do;
  eq[36] := simplify(add(w5 [i]*b[i]*c[i], i=3..s));
  #print(`  eq[36]====`, eq[36]);
  
  Dnx := 4*c[3]^2*(-4*c[4]^2+5*c[3]*c[4]+3*c[5]*c[4]-3*c[3]*c[5])*a[5,2]
        -3*c[5]^2*(3*c[3]-2*c[4])*(6*c[3]*c[4]-4*c[3]*c[5]-4*c[5]*c[4]+3*c[5]^2);
  x4  := 3/4*((-4*c[3]^2*(c[3]*c[4]+3*c[3]*c[5]-2*c[4]^2))*a[5,2]
        +3*c[5]^2*(3*c[3]-2*c[4])*(3*c[3]*c[4]+c[3]*c[5]-2*c[5]*c[4]))/c[3]^2/Dnx;
  x2  := -4/9*c[3]^2*x4;
  x4a := 4*c[3]^2*c[4]*(-c[4]+2*c[3])*a[5,2] -3*c[5]^2
          *(9*c[3]^2*c[4]-3*c[3]^2*c[5]-3*c[4]^2*c[3]-7*c[3]*c[5]*c[4]+3*c[5]^2*c[3]+2*c[5]*c[4]^2);
  x4a := x4a/Dnx;
  x3  := 1/27*c[3]*(4*x4*c[3]^2-27*x4a);

  # Z=a62*c2^3 + a63*c3^3 + a64*c4^3 + a65*c5^3
  # t:=simplify(x2*w2+x3*w31+x4*w4+x4a*Aew31-Aw31): 777;

  Zt:=a[6,2]*c[2]^3 + a[6,3]*c[3]^3 + a[6,4]*c[4]^3 + a[6,5]*c[5]^3;
  #eq[37]:= x3*a[6,2] + x4*(4*Zt-c[6]^4) + x4a*a[6,2]*c[6]-add(a[6,i]*a[i,2], i=3..5);
  eq[38] := x3*a[7,2] + x4*w4s7          + x4a*a[7,2]*c[7]-add(a[7,i]*a[i,2], i=3..6);
  eq[39] := x3*a[8,2] + x4*w4s8          + x4a*a[8,2]*c[8]-add(a[8,i]*a[i,2], i=3..7);
  eq[40] := x3*a[9,2] + x4*w4s9          + x4a*a[9,2]     -add(a[9,i]*a[i,2], i=3..8);

  #-------------------------------------------------------------------------------------------------

  sol:=solve({eq[21],eq[28],eq[29],eq[38]},{a[7,3], a[7,4],a[7,5], a[7,6]});
  assign(sol);

  sol:=solve({eq[22],eq[23],eq[30],eq[31],eq[32], eq[33],eq[36],eq[39],eq[40]},
             {a[8,3],a[8,4],a[8,5],a[8,6],a[8,7], a[9,3],a[9,4],a[9,5],a[9,6]});
  assign(sol);


  # ---------- finish --------------------------------------------- 
  A := Matrix(10,10);
  for i from 2 to s do for j from 1 to i-1 do A[i,j]:=simplify(eval(a[i,j])); end do; end do; 
  for j from 1 to s do A[s+1,j]:=simplify(eval(b[j])); end do;
  for i from 2 to s+1 do
    A[i,1]:=simplify(eval(c[i] - add(A[i,j], j=2..i-1)));
  end do;
  simplify(A);

end proc:

# ==========================================================
RK_79b := proc(c2,c3,c4,c5,a52,c7,c8)::Matrix; local A, n, c, i, k, lam, s, an;
  A := RK_79a(c3,c4,c5,a52,c7,c8);
  n :=RowDimension(A);
  lam := c2/A[2,1];
  A[2,1]:=c2;
  for k from 3 to n do
      s := A[k,1]+A[k,2];
      A[k,2]:=A[k,2]/lam;
      A[k,1] := s - A[k,2]; 
  end do;


  simplify(A);
end proc:
# ==========================================================
RK_79c := proc(c3,c4,c5,a52,c7, ind)::Matrix; local c8, eq78, A;
  eq78:=RK79_c7c8(c3,c4,c5,a52);
  #print(777, eq78);
  c8 :=solve(eq78(c7,c8), c8)[ind];
  #print(778, c8);
  A := RK_79a(c3,c4,c5,a52,c7,c8);
  simplify(A);
end proc:
# ==========================================================
RK_79d := proc(c2, c3,c4,c5,a52,c7, ind)::Matrix; local c8, eq78, A;
  eq78:=RK79_c7c8(c3,c4,c5,a52);
  c8 :=solve(eq78(c7,c8), c8)[ind];
  A := RK_79b(c2,c3,c4,c5,a52,c7,c8);
  simplify(A);
end proc: