BeginPackage["RungeKuttaNDSolve`"]
RungeKuttaNDSolve::usage = "RungeKuttaNDSolve[deqn,vars,{t,t0,t1}, opts] \
solve a first order autonomous system of ordinary differential equations \
deqn, depending on the list of variables vars on the interval t in [t0,t1]. \
The dependence on t must be given explicit in the variables. The initial \
values for the variables at t==t0 must be included int the list of equations \
deqn. The RKReturn option determines the possible results. RKReturn -> \
RKEndPoint return the value of the vars at t1 (and not an interpolation). \
RKReturn->RKInterpolation returns an interpolation function like NDSolve[] \
does. RKReturn->RKPoincareSection return the points where the solution cross \
the hyperplane given by the PSHyperplane option. The RKStoppingFunction \
option can be used to terminate the integration before t1 is reached. \
StartingStepSize, MaxSteps, PrecisionGoal and MaxStepSize options can be used \
in the same way as in NDSolve[]."

rk43Interpolation::usage = "rk43Interpolation[u,y,h] is the continuous \
interpolation for the RKInterpolateLow method. It interpolates y[t+u*h] in \
the current interval. It should only called in the function given in the \
RKOutputFunction option to RungeKuttaNDSolve[]. "

rk54DPInterpolation::usage = "rk54DPInterpolation[u,y,h] is the continuous \
interpolation for the RKInterpolationMedium method.It interpolates y[t+u*h] \
in the current interval. It should only called in the function given in the \
RKOutputFunction option to RungeKuttaNDSolve[]."

rk65Interpolation::usage = "rk65Interpolation[u,y,h] is the continuous \
interpolation for the RKInterpolation method.It interpolates y[t+u*h] in the \
current interval. It should only called in the function given in the \
RKOutputFunction option to RungeKuttaNDSolve[]."

rk87DPInterpolation::usage = "rk87DPInterpolation[u,y,h] is the continuous \
interpolation for the RKInterpolationHigh method.It interpolates y[t+u*h] in \
the current interval. It should only called in the function given in the \
RKOutputFunction option to RungeKuttaNDSolve[]."

(* Methods used by RungeKuttaNDSolve[] *)

RKInterpolate::usage = "RKInterpolate is a possible value of the Method \
option of RungeKuttaNDSolve[], is uses a 6th order Verner method with \
continuous output."

RKInterpolateLow::usage = "RKInterpolateLow is a possible value of the Method \
option of RungeKuttaNDSolve[], is uses a 3th order Verner method with \
continuous output."

RKInterpolateMedium::usage = "RKInterpolateMedium is a possible value of the \
Method option of RungeKuttaNDSolve[], is uses a 5th order Dormand/Prince \
method with continuous output."

RKStability::usage = "RKStability is a possible value of the Method option of \
RungeKuttaNDSolve[], is uses a 5th order method
with a large stability interval."

RKInterpolateHigh::usage = "RKInterpolateHigh is a possible value of the \
Method option of RungeKuttaNDSolve[], is uses a 8th order Dormand/Prince \
method with continuous output."

RKPhaseLag::usage = "RKPhaseLag is a possible value of the Method option of \
RungeKuttaNDSolve[], is uses a 5th order method with a small error in the \
phase for oscillating systems."

(* Additional options for RungeKuttaNDSolve[] *) 

(#::"usage" = 
        ToString[#] <> 
          "is an option to RungeKuttaNDSolve[]. It must be a function with \
the parameters t, y[t], f[y],x+h, y[t+h], f[t[t+h], h in this order. It us \
called after every successful step. It should return True or Null to the \
solver to continue the integration, otherwise the solver stops and returns \
the result.") & /@ {RKOutputFunction, RKStoppingFunction}

RKReturn::usage = "RKReturn is an option to RungeKuttaNDSolve[]. It \
determines the returned result, the option value can be RKEndPoint (default), \
RKInterpolation or RKPoincareSection."
  
RKStatistics::usage = "RKStatistics is an option to RungeKuttaNDSolve[]. It \
prints an Message[] about the number of steps, the number of evaluations of \
the right hand side and the number of fail steps."

RKInterpolationOrder::usage = "RKInterpolationOrder is an option to \
RungeKuttaNDSolve[]. It should be an integer or Automatic. By default \
it creates an interpolation of the order of the method."

PSHyperplane::usage = "PSHyperplane is an option to RungeKuttaNDSolve[]. It \
should be an equation that must be valid for the solution points. When the \
RKReturn->RKPoincareSection option is given, RungeKuttaNDSolve[] look for \
points along the solution of the system where the equation is valid."

RKEndPoint::usage = "RKEndPoint is a possible value for the RKReturn option \
of RungeKuttaNDSolve[]."

RKInterpolation::usage = "RKInterpolation  is a possible value for the RKReturn option of RungeKuttaNDSolve[]."
RKPoincareSection::usage = "RKPoincareSection  is a possible value for the RKReturn option of RungeKuttaNDSolve[]."

Begin["`Private`"]

(*
  The following two methods are from
    "A Family of Fifth-Order Runge-Kutta Pairs"
    S. N. Papakostas and G. Papageorgiou
    Mathematics of Computation, Vol. 65, No 215, July 1996, 
  Pages 1165 - 1181 
  *)
  
normn = Compile[{{x, _Real, 1}}, Sqrt[Dot[x, x]]]

rkStepPhase56[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
        b21 = 21/83,
        b31 = 1208/11109, b32 = 2656/11109,
        b41 = 18692077/12446784, b42 = -12162239/1555848, 
        b43 = 4392287/592704,
        b51 = 208414811/159752250, b52 = -4726684/694575, 
        b53 = 783101447/119400750, b54 = 8928/1037875,
        b61 = 22441327/25440576, b62 = -631879/138264, 
        b63 = 2740042785/589452352, b64 = -652680/6833369, b65 = 55125/408208,          
        c = h*{919/8832, 17070301/32230080, 180075/66424, -18375/3968, 823/360},   
        cc = h*{9753/106720, 148473901/259631200, 3858407/9631480, -15925/14384, 5761/5800, 1/20}},
      Block[{delta, err, hfac, y1},
        $k1 = f1;
        $k2 = f[y + b21*h*$k1];
        $k3 = f[y + h*(b31*$k1 + b32*$k2)];
        $k4 = f[y + h*(b41*$k1 + b42*$k2 + b43*$k3)];
        $k5 = f[y + h*(b51*$k1 + b52*$k2 + b53*$k3 + b54*$k4)];
        (* instead of creating a new $k6 we overwrite $k2 *)
        $k2 = f[y + h*(b61*$k1 + b62*$k2 + b63*$k3 + b64*$k4 + b65*$k5)];
        y1 = y + Dot[c, {$k1, $k3, $k4, $k5, $k2}];
        $k7 = f[y1];
        $FunctionCalls += 6;
        delta = y1 - y - Dot[cc, {$k1, $k3, $k4, $k5, $k2, $k7}];
        err = Max[normn[delta], $MachineEpsilon];
        hfac = saveFac*
          If[err < locError, (locError/err)^(1/6), (locError/err)^(1/5)];
        {y1, $k7, Min[Max[1/10, hfac], 5], err}
      ]
    ]

rkStepMaxStability56[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
        b21 = 6/29,
        b31 = 225/3136, b32 = 783/3136,
        b41 = 239628760/2122043913, b42 = 25435900/707347971, 
        b43 = 684353600/2122043913,
        b51 = 1280931436/7518984375, b52 = -187500341/916600000, 
        b53 = 4930503382/9237609375, b54 = 485982986841/6897415000000,
        b61 = -1647119/2711961, b62 = 1773611/1549692, 
        b63 = 49937925758/14010377949, b64 = -4391450908515/489152454812, 
        b65 = 193373956250/32908097043,
        c = h*{881/12312, 54762008/55654857, -7780827681/4587493528, 2864375000/1873713519, 14349/128269}, 
        cc = h*{386231063/3718839600, 153892305602/240150707955, -1301184381146067/1385652420132400, 1721174293750/1617014766897, 44252316/553480735, 1/20}},
      Block[{delta, err, hfac, y1},
        $k1 = f1;
        $k2 = f[y + b21*h*$k1];
        $k3 = f[y + h*(b31*$k1 + b32*$k2)];
        $k4 = f[y + h*(b41*$k1 + b42*$k2 + b43*$k3)];
        $k5 = f[y + h*(b51*$k1 + b52*$k2 + b53*$k3 + b54*$k4)];
        (* instead of creatng a new $k6 we overwrite $k2 *)
        $k2 = f[y + h*(b61*$k1 + b62*$k2 + b63*$k3 + b64*$k4 + b65*$k5)];
        y1 = y + Dot[c, {$k1, $k3, $k4, $k5, $k2}];
        $k7 = f[y1];
        $FunctionCalls += 6;
        delta = y1 - y - Dot[cc, {$k1, $k3, $k4, $k5, $k2, $k7}];
        err = Max[normn[delta], $MachineEpsilon];
        hfac = saveFac*
          If[err < locError, (locError/err)^(1/6), (locError/err)^(1/5)];
        {y1, $k7, Min[Max[1/10, hfac], 5], err}
      ]
    ]  


(* y  - y[n]  solution at t[n] *)
(* f1 - f[y[n]]  the right hand side of the ode´s at y[n] *)
(* f  - the Function[], CompiledFunction[] for the right hand sides *)
(* h  - the integration step *)
(* locErr - the desired local error *) 
rkStep34[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
        b21 = 3/8,
        b31 = 9/64, b32 = 27/64,
        b41 = 925/5184, b42 = 625/3456, b43 = 4375/10368,
        b51 = 63863/135675, b52 = -2516/1809, b53 = 11872/5427, b54 = -448/1675,
        c = h*{9./50, 128/147, -1024/3675, 67/294},
        d = h*{-67/1350, 536/1323, -2144/3675, 67/294}},
    Block[{delta, err, hfac, y1},
      $k1 = f1;
      $k2 = f[y + b21*h*$k1];
      $k3 = f[y + h*(b31*$k1 + b32*$k2)];
      $k4 = f[y + h*(b41*$k1 + b42*$k2 + b43*$k3)];
      $k2 = f[y + h*(b51*$k1 + b52*$k2 + b53*$k3 + b54*$k4)];
      y1 = y + Dot[c, {$k1, $k3, $k4, $k2}];
      $FunctionCalls += 5;
      delta = Dot[d, {$k1, $k3, $k4, $k2}];
      err = Max[normn[delta], $MachineEpsilon];
      hfac = 
        saveFac*If[err < locError, (locError/err)^(1/4), (locError/err)^(1/3)];
      {y1, f[y1], Min[Max[1/10, hfac], 5], err}
     ]
   ]

RungeKuttaNDSolve::nointp = "Call to continuous output function must be done during a Runge-Kutta step."

(* Global values $RKActive, $k1, $k2, $k3, $k4 *)
rk34Interpolation[u_, y_, h_] :=
  Module[{u2 = u^2, b1, b3, b4, b5},
    If[! $RKActiveQ,
      Message[ RungeKuttaNDSolve::nointp];
      Return[y]
     ];
    b1 = h*u*(1 + u*(-913/450 + (16/9 - (128*u)/225)*u));
    b3 = h*(3200/441 + u*(-4864/441 + (2048*u)/441))*u2;
    b4 = h*(-9216/1225 + (2048/147 - (8192*u)/1225)*u)*u2;
    b5 = h*(225/98 + u*(-688/147 + (128*u)/49))*u2;
    y + (b1*$k1 + b3*$k3 + b4*$k4 + b5*$k2)
   ]
   
rkStepVerner65[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
    b21 = 1/8,
    b31 = 1/18, b32 = 1/9,
    b41 = 1/16, b42 = 0, b43 = 3/16,
    b51 = 1/4, b52 = 0, b53 = -3/4, b54 = 1,
    b61 = 134/625, b62 = 0, b63 = -333/625, b64 = 476/625, b65 = 98/625,
    b71 = -98/1875, b72 = 0, b73 = 12/625, b74 = 10736/13125, b75 = -1936/1875, b76 = 22/21,
    b81 = 9/50, b82 = 0, b83 = 21/25, b84 = -2924/1925, b85 = 74/25, b86 = -15/7, b87 = 15/22,
    c = h*{11/144, 256/693, 125/504, 125/528, 5/72},
    cc = h*{1/18, 32/63, -2/3, 125/126, -5/63, 4/21}     
   },
   Block[{delta, err, hfac, y1},
     $k1 = f1;
     $k2 = f[y + b21*h*$k1];
     $k3 = f[y + h*(b31*$k1 + b32*$k2)];
     $k4 = f[y + h*(b41*$k1 + b43*$k3)];
     $k5 = f[y + h*(b51*$k1 + b53*$k3 + $k4)];
     $k6 = f[y + h*(b61*$k1 + b63*$k3 + b64*$k4 + b65*$k5)];
     $k7 = f[y + h*(b71*$k1 + b73*$k3 + b74*$k4 + b75*$k5 + b76*$k6)];
     $k8 = f[y + h*(b81*$k1 + b83*$k3 + b84*$k4 + b85*$k5 + b86*$k6 + b87*$k7)]; 
     y1 = y + Dot[c, {$k1, $k4, $k6, $k7, $k8}];
     $k9 = f[y1];
     $FunctionCalls += 8;
     delta = y1 - y - Dot[cc, {$k1, $k4, $k5, $k6, $k8, $k9}];
     err = Max[normn[delta], $MachineEpsilon];
     hfac = saveFac*
          If[err < locError, (locError/err)^(1/7), (locError/err)^(1/6)];
     $RKMoreKiPresentQ = False;
     {y1, $k9, Min[Max[1/10, hfac], 5], err}
    ]
  ]

rk65Interpolation[u_, y_, h_] :=
  ( If[!$RKActiveQ,
      Message[ RungeKuttaNDSolve::nointp];
      Return[y]
     ];
    If[! $RKMoreKiPresentQ,(* The compute the additional k's *)
       With[{
         b10 = h*{173/2304, 8/21, -1/9, 1625/8064, -125/2304, -445/8064, 127/2016},
         b11 = h*{153/2048, 237/616, -33/128, 375/896, -375/22528, -45/896, 201/3584, 9/64},
         b12 = h*{60201/800000, 91368/240625, -81/625, 8073/22400, 10449/70400, 13851/560000, -7533/700000, 162/3125}},
         $k10 = $RKRighthandSide[y + Dot[b10, {$k1, $k4, $k5, $k6, $k7, $k8, $k9}]];
         $k11 = $RKRighthandSide[y + Dot[b11, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10}]];
         $k12 = $RKRighthandSide[y + Dot[b12, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10}]];
         $FunctionCalls += 3;
        ];
       $RKMoreKiPresentQ = True;
      ]; 
    With[{bri = {
      {-2125/1053, 5771/702, -37577/2808, 92843/8424, -8921/1872, 1}, 
      {102400/9009, -6144/143, 186880/3003, -17920/429, 11520/1001, 0}, 
      {6250/819, -375/13, 45625/1092, -4375/156, 5625/728, 0}, 
      {3125/429, -7875/286, 45625/1144, -30625/1144, 16875/2288, 0}, 
      {250/117, -105/13, 1825/156, -1225/156, 225/104, 0}, 
      {200/13, -548/13, 536/13, -215/13, 27/13, 0}, 
      {-400/13, 1408/13, -3691/26, 1067/13, -459/26, 0}, 
      {12800/351, -12800/117, 13376/117, -16256/351, 64/13, 0}, 
      {-50000/1053, 50000/351, -109375/702, 78125/1053, -3125/234, 0}}},
       y + h*Dot[
            u*(Fold[(u*#1 + #2) &, First[#], Rest[#]] & /@ bri), 
            {$k1, $k4, $k6, $k7, $k8, $k9, $k10, $k11, $k12}]
     ])   

(* 
  Dormand, Prince methods. Numerical values of the coefficients are from
      ftp : // ftp.tees.ac.uk/pub/j.r.dormand/rkcoeff.dat
  *)
  
rkStepDormandPrince54[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
    b2 = h/5,
    b3 = h*{3/40, 9/40},
    b4 = h*{44/45, -56/15, 32/9},
    b5 = h*{19372/6561, -25360/2187, 64448/6561, -212/729},
    b6 = h*{9017/3168, -355/33, 46732/5247, 49/176, -5103/18656},
    b7 = h*{35/384, 0, 500/1113, 125/192, -2187/6784, 11/84},
    c = h*{35/384,  500/1113, 125/192, -2187/6784, 11/84},
    cc = h*{5179/57600, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40}},
    Block[{delta, err, hfac, y1},
      $k1 = f1;
      $k2 = f[y + b2*$k1];
      $k3 = f[y + Dot[b3, {$k1, $k2}]];
      $k4 = f[y + Dot[b4, {$k1, $k2, $k3}]];
      $k5 = f[y + Dot[b5, {$k1, $k2, $k3, $k4}]];
      $k6 = f[y + Dot[b6, {$k1, $k2, $k3, $k4, $k5}]];
      y1 = y + Dot[c, {$k1, $k3, $k4, $k5, $k6}];
      $k7 = f[y1];
      $FunctionCalls += 6;
      delta = y1 - y - Dot[cc, {$k1, $k3, $k4, $k5, $k6, $k7}];
      err = Max[normn[delta], $MachineEpsilon];
      hfac = saveFac*
          If[err < locError, (locError/err)^(1/6), (locError/err)^(1/5)];
      $RKMoreKiPresentQ = False;
      {y1, $k7, Min[Max[1/10, hfac], 5], err}
     ]
   ]

rk54DPInterpolation[u_, y_, h_] :=
  ( If[!$RKActiveQ,
      Message[ RungeKuttaNDSolve::nointp];
      Return[y]
     ];
    If[! $RKMoreKiPresentQ,(* The compute the additional k's *)
      With[{
        b8 = h*{6245/62208, 8875/103032, -125/1728, 801/13568, -13519/368064, 11105/368064},              
        b9 = h*{632855/4478976,  4146875/6491016, 5490625/14183424, -15975/108544, 8295925/220286304, -1779595/62938944, -805/4104}},
        $k8 = $RKRighthandSide[y + Dot[b8, {$k1, $k3, $k4, $k5, $k6, $k7}]];
        $k9 = $RKRighthandSide[y + Dot[b9, {$k1, $k3, $k4, $k5, $k6, $k7, $k8}]];
        $FunctionCalls += 2;
        ];
      $RKMoreKiPresentQ = True;
     ]; 
    With[{bri = {{3303/880, -19683/1760, 125923/10560, -38039/7040, 1},
                 {36000/4081, -90000/4081, 205000/12243, -12500/4081, 0},
                 {1125/88, -5625/176, 25625/1056, -3125/704,0},
                 {-59049/9328, 295245/18656, -448335/37312, 164025/74624, 0},
                 {18/7, -45/7, 205/42, -25/28, 0},
                 {108/55, -171/55, 73/55, -2/11, 0},
                 {-648/55, 3537/110, -1593/55, 189/22, 0},
                 {-648/55, 2943/110, -999/55, 351/110, 0}}},
       y + h*Dot[
               u*(Fold[(u*#1 + #2) &, First[#], Rest[#]] & /@ bri), 
               {$k1, $k3, $k4, $k5, $k6, $k7, $k8, $k9}
              ]
     ])

rkStepDormandPrince87[y_, f1_, f_, h_, locError_] :=
  With[{saveFac = 9/10,
    b2 = h*832108848095584489/15819514147941225301, 
    b3 = h*{2496326544286753467/126556113183529802408, 7488979632860260401/126556113183529802408},     
    b4 = h*{841554882220704499/28442837381251792786, 2524664646662113497/28442837381251792786},
    b5 = h*{2861105270853988/11853846583185727, -1121928579081792012/1268361584400872789, 1173023921684928458/1268361584400872789},
    b6 = h*{1/27, 194058014822681543/1135980771991148295,142528880582103137/1135980771991148295},
    b7 = h*{19/512,  443991049/2607843072, 39258821/651960768, -9/512},
    b8 = h*{13772/371293,  417418705735041635/2449871394793524377, 
            20213706164464769/188451645753348029, -5688/371293, 3072/371293},
    b9 = h*{58656157643/93983540625, -24785085474788602622617/
              7374554384043849721875, -32013653952816691272967/
              36872771920219248609375, 96044563816/3480871875, 
            5682451879168/281950621875, -165125654/3796875},
    b10 = h*{8909899/18653125, -325643933915917646309371/
             130879796011819495500000, -77257143004054276420613/
             130879796011819495500000, 96663078/4553125, 
              2107245056/137915625, -4913652016/147609375, -78894270/3880452869},
    b11 = h*{-20401265806/21769653311, 
            38171323495890517387914273/7359927197591963973146245, 
            8032899429348702276090447/7359927197591963973146245, -43306765128/
              5313852383, -20866708358144/1126708119789, 
            14886003438020/654632330667, 
            35290686222309375/14152473387134411, -1477884375/485066827},
    b12 = h*{39815761/17514443, -489691433455557804808299373/
              46484564448460036998778525, -93009666143676504491807127/
              46484564448460036998778525, -844554132/47026969, 
            8444996352/302158619, -2509602342/877790785, -28388795297996250/
              3199510091356783, 226716250/18341897, 1371316744/2131383595},
    cc = h*{1766245085361/27894533582822, 2, 
            31667627311137/27894533582822, -37871782384591/13947266791411, 
            1506508387605/27894533582822, 2885889429784/13947266791411, 
            2944844885381/13947266791411, 1247181771927/27894533582822},
    c = h*{104257/1920240, 3399327/763840, 
            66578432/35198415, -1674902723/288716400, 
            54980371265625/176692375811392, -734375/4826304, 
            171414593/851261400, 137909/3084480}},
    Block[{delta, err, hfac, y1},
      $k1 = f1;
      $k2 = f[y + b2*$k1];
      $k3 = f[y + Dot[b3, {$k1, $k2}]];
      $k4 = f[y + Dot[b4, {$k1, $k3}]];
      $k5 = f[y + Dot[b5, {$k1, $k3, $k4}]];
      $k6 = f[y + Dot[b6, {$k1, $k4, $k5}]];
      $k7 = f[y + Dot[b7, {$k1, $k4, $k5, $k6}]];
      $k8 = f[y + Dot[b8, {$k1, $k4, $k5, $k6, $k7}]];
      $k9 = f[y + Dot[b9, {$k1, $k4, $k5, $k6, $k7, $k8}]];
      $k10 = f[y + Dot[b10, {$k1, $k4, $k5, $k6, $k7, $k8, $k9}]];
      $k11 = f[y + Dot[b11, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10}]];
      $k12 = f[y + Dot[b12, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10, $k11}]];
      y1 = y + Dot[c, {$k1, $k6, $k7, $k8, $k9, $k10, $k11, $k12}];  
      delta = y1 - y - Dot[cc, {$k1, $k6, $k7, $k8, $k9, $k10, $k11, $k12}];
      $FunctionCalls += 12;
      err = Max[normn[delta], $MachineEpsilon]; 
      hfac = saveFac*
      If[err < locError, (locError/err)^(1/8), (locError/err)^(1/7)];
      $RKMoreKiPresentQ = False;
      {y1, f[y1], Min[Max[1/10, hfac], 5], err}
     ]
   ]

rk87DPInterpolation[u_, y_, h_] :=
  ( If[! $RKActiveQ,
      Message[ RungeKuttaNDSolve::nointp];
      Return[y]
     ];
    If[! $RKMoreKiPresentQ,(* The compute the additional k's *)
      With[{
        b13 = h*{104257/1920240, 3399327/763840, 
                66578432/35198415, -1674902723/288716400, 
                54980371265625/176692375811392, -734375/4826304, 
                171414593/851261400, 137909/3084480},
        b14 = h*{61937212127458037412818875/2173408747461889498729075762, 
                488968626810081815037897585/2173408747461889498729075762, 
                42593854112307271975723718/
                  1086704373730944749364537881, -1283656793492287688396290190/
                  1086704373730944749364537881, -780230339235921612430671017/
                  1086704373730944749364537881, 
                3784213786917400917091588683/
                  2173408747461889498729075762, -76512361440119541449704287/
                  2173408747461889498729075762, 
                63161319624717787862989590/1086704373730944749364537881, 
                162258289927658535869076/
                  1086704373730944749364537881, -2841389793059711294432392/
                  1086704373730944749364537881, 
                3092359402318008547543084/1086704373730944749364537881},
        b15 = h*{-14185244866816877137470587/1815617987128373777220360800, 
                2949944204646799624765576873/5446853961385121331661082400, 
                47424215805529698253817111/
                  340428372586570083228817650, -4788590462337322066667487151/
                  1361713490346280332915270600, -3542138193965218116021530583/
                  1815617987128373777220360800, 
                2318928771095344872636955577/
                  453904496782093444305090200, -522098424964402205454951287/
                  5446853961385121331661082400, 
                219720821463750588913510553/1361713490346280332915270600, 
                116669178953177948544727/
                  340428372586570083228817650, -3710177330343145593647757/
                  453904496782093444305090200, 
                8005663818985003345755107/
                  907808993564186888610180400, -792427964128249312487849759/
                  5446853961385121331661082400},
       b16 = h*{-253246664612804731174409126/
                  4298318608082284612088123043, 
                412378479461865733055566529/452454590324451011798749794, 
                3751593077197114123551216349/
                  8596637216164569224176246086,-48259851414567510834337979773/
                  8596637216164569224176246086, -799567868155999074408319364/
                  204681838480108791051815383, 73704621495236313126185799035/
                  8596637216164569224176246086, -203812829974244903044089515/
                  1228091030880652746310892298, 1152718146300914113953218405/
                  4298318608082284612088123043, 7941461350160774510596889/
                  1432772869360761537362707681, -1501234390372394975966639/
                  330639892929406508622163311, 6203701716626364228473688/
                  1432772869360761537362707681, -95312072573093119921003097/
                  330639892929406508622163311, 2/7}
      },
     $k13 = $RKRighthandSide[y + Dot[b13, {$k1, $k6, $k7, $k8, $k9, $k10, $k11, $k12}]];
     $k14 = $RKRighthandSide[y + Dot[b14, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10, $k11, $k12, $k13}]];
     $k2 =  $RKRighthandSide[y + Dot[b15, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10, $k11, $k12, $k13, $k14}]];
     $k3 =  $RKRighthandSide[y + Dot[b16, {$k1, $k4, $k5, $k6, $k7, $k8, $k9, $k10, $k11, $k12, $k13, $k14, $k2}]];
     $FunctionCalls += 4;
    ];
   $RKMoreKiPresentQ = True;
  ]; 
  With[{bri = {{31209862049273161473565530870/1188057472418147903718890477, 
               -398811383198491022942114013049/4196488651810448704008105265, 
              1104720771870777038977309621664/8048101236095571817892380181, 
              -295916641860596250213324720927/2908501044317972075231454589, 
              3842358082418320312630980939/92875772330122244066807290, 
              -34955677982384479519902809081/3852053886053326660818139193, 1}, 
              {-173795253237022764697084601029/1188057472418147903718890477, 
              1682173749941018145613355622063/4196488651810448704008105265, 
              -1740760865239756105143887727578/8048101236095571817892380181, 
              -699987094255637099202075935138/2908501044317972075231454589, 
              138944356898435219858768150509/557254633980733464400843740, 
              -163675568899079121484659202504/3852053886053326660818139193, 0}, 
              {1342550138048327712351704695565/1188057472418147903718890477, 
               -3210975227783050730204109002001/839297730362089740801621053, 
              40266183855064513864648636375875/8048101236095571817892380181,
                -9191716014386437947320800339405/2908501044317972075231454589, 
              267536207629801821165117790187/278627316990366732200421870, 
              -406197406808662645735088627563/3852053886053326660818139193, 0}, 
              {-835376209896779338701047834774/1188057472418147903718890477, 
              10095374043740569932173641308249/4196488651810448704008105265, 
              -27016870126160508602294548475275/8048101236095571817892380181, 
              7106851197722504370079887793960/2908501044317972075231454589, 
              -126989732840120246719331849621/139313658495183366100210935, 
              449377370612179412357263463031/3852053886053326660818139193, 0}, 
              {118523914171977715788813075462/1188057472418147903718890477, 
              -1440500437284736806332951089017/4196488651810448704008105265, 
              3587401253682266076969938128314/8048101236095571817892380181, 
              -779954882531904914239657230116/2908501044317972075231454589, 
              3376474218013923070601985298/46437886165061122033403645, 
              -2268857000812900165739251218/350186716913938787347103563, 0}, 
              {-137721017483539677801021223384/1188057472418147903718890477, 
              1723552682402200084433680927401/4196488651810448704008105265, 
                -4442757322820882742957783954256/8048101236095571817892380181, 
              1002934995236068121159515604755/2908501044317972075231454589, 
              -53787945977582342180299707431/557254633980733464400843740, 
              33804704888800897793969090552/3852053886053326660818139193, 0}, 
              {-27098195885822268859586193958/1188057472418147903718890477, 
              274982070053409786826517201537/4196488651810448704008105265, 
              -566303954209189436198588528212/8048101236095571817892380181, 
              103197190772348204543369503390/2908501044317972075231454589, 
              -4643676787908105868482084017/557254633980733464400843740, 
              2697614073343118802856591919/3852053886053326660818139193, 0}, 
              {-10797824299363167239600282630/1188057472418147903718890477, 
              107216461668477853975060899906/4196488651810448704008105265, 
              -210889497375033136609278678801/8048101236095571817892380181, 
              35327576926576090982086258626/2908501044317972075231454589, 
              -468896405834309472421535429/185751544660244488133614580, 
              639041948989764734497745437/3852053886053326660818139193, 0}, 
              {17902029353873180247342195331/1188057472418147903718890477, 
              -173109236120268042243256758081/4196488651810448704008105265, 
              336989913597520940480710867953/8048101236095571817892380181, 
              -56910811838915522436594842075/2908501044317972075231454589, 
              2329590455306601148108390913/557254633980733464400843740,
              -1166182759577769747292967340/3852053886053326660818139193, 0}, 
              {-369052421188973545981266240922/1188057472418147903718890477, 
              4479361905873801469016720817228/4196488651810448704008105265,
              -11524687062474207127893844010983/8048101236095571817892380181, 
              2753790668520802284568893959329/2908501044317972075231454589, 
              -29175999946430909363399346829/92875772330122244066807290, 
              163852753845995734458830288048/3852053886053326660818139193, 0}, 
              {-187012177416630123939634651424/1188057472418147903718890477, 
              2331527890304218032226862253589/4196488651810448704008105265, 
              -5860920228333332676872581154805/8048101236095571817892380181, 
              1231177150560210276971221974242/2908501044317972075231454589, 
              -13542301973509461899033528567/139313658495183366100210935, 
              15256292127344811892242896721/3852053886053326660818139193, 0}, 
              {230667155784679117357815530893/1188057472418147903718890477, 
              -2626891608464945781726972159821/4196488651810448704008105265, 
              6067893262397831906893917536104/8048101236095571817892380181,
              -1208793334865017614892522026641/2908501044317972075231454589, 
              11420849838415041011956432091/111450926796146692880168748, 
              -34675514038007821729584705822/3852053886053326660818139193, 0}}},
       y + h*Dot[
               u*(Fold[(u*#1 + #2) &, First[#], Rest[#]] & /@  bri), 
               {$k1, $k6, $k7, $k8, $k9, $k10, $k11, $k12, $k13, $k14, $k2, $k3}
              ]
      ])
      
$ContinueSym = {Null, True};

rkDSolve[y0_, func_, hInit_?NumericQ, {x1_?NumericQ, x2_?NumericQ},
     localErr_, maxsteps_Integer, t_, rkstep_, output_:None, stop_:None] :=
  Module[{y = y0, f, x = x1, h, y1, f1, errorEst, hfac, hnew, succ, res, i = 0},
    If[(x2 - x1)*h < 0, Return[$Failed]];
    hnew = hInit;
    $RKRighthandSide = func; (* for the additional function evaluation in the high order Verner65 method *) 
    f = func[y];
    If[None =!= output, output[x, y, f, x, y, f, 0]]; (* Call the output function for the initial values *) 
    If[None =!= stop, stop[x, y, f, x, y, f, 0]];
    $RKActiveQ = True;
    While[Abs[x2 - x] > $MachineEpsilon,
      h = hnew;
      If[Abs[h] < $MachineEpsilon, (* Step size too small ... *)
        Message[NDSolve::ndsz, t, x];
        Break[];
       ];
      $TotalSteps++;
      If[Abs[x + h - x1] > Abs[x2 - x1], h = x2 - x]; (* Don't overshoot the end *)
      (* Integrate one step *)
      {y1, f1, hfac, errorEst} = rkstep[y, f, func, h, localErr];
      hnew = hfac*h;
      If[TrueQ[errorEst < localErr],(* Accept the step, increment the variables *)
        If[None =!= output  && !MemberQ[$ContinueSym, res = output[x, y, f, x + h, y1, f1, h]],
          Break[];
         ];  
        If[None =!= stop && !MemberQ[$ContinueSym, res = stop[x, y, f, x + h, y1, f1, h]],
          Break[];
         ];
        x += h;
        y = y1;
        f = f1, 
      (* Else *)
        $FailSteps++
       ];     
      If[++i >= maxsteps,
        Message[NDSolve::mxst, maxsteps, t, x];
        Break[]
       ];
     ]; (* While *)
    $RKActiveQ = False;
    $RKMoreKiPresentQ = False; 
    y1
    ] (* Module *)      
  
  
Clear[makeODEFunction]

makeODEFunction[f_List, vars_, compile_] :=
  Module[{ff, sym, newsym, rules},
    sym = Unique[];
    ff = f /. Thread[
                vars -> (Hold[Part[newsym, #]] & /@ Table[i, {i, Length[vars]}] /.
                            newsym -> sym)
               ];
    If[compile,
       Hold[Evaluate[ {{{sym, _Real, 1}}, ff}] ] /.  
          Hold[a_Part] :> a /. Hold[l_] :> Apply[Compile, Unevaluated[l]],
    (* Else *)
      ff = Hold[Evaluate[ff]] /. Hold[a_Part] :> a;
      ReleaseHold[Function[Evaluate[{sym}], #] & @@ {ff}]
     ]
   ]
   
FirstOrderSystemQ[deqn_, vars_, t_] :=
    And @@ (#1 == 1 & /@ (Flatten[
                             Cases[
                               deqn, 
                               Derivative[_][Head[#]][t], 
                               Infinity] & /@ vars
                            ] /. 
                               Derivative[n_][_][t] :> n))
                               
AutonomDEqnQ[deqn_, vars_List, t_] :=
  Module[{rep},
    rep = 
      (# -> (# /. a_[t] :> a) ) & /@ (Join @@ ({#, D[#, t]} & /@ vars));
    FreeQ[deqn /. rep, t]
   ]                               

AutonomDEqnQ[deqn_, vars_, t_] :=
  Module[{rep},
    rep = 
      (# -> (# /. a_[t] :> a) ) & /@ ({#, D[#, t]} & [vars]);
    Print[vars,"  ",rep];
    FreeQ[deqn /. rep, t]
   ]  

makeODESystem[sys_, vars_, t_, t0_, compile_:True] :=
  Module[{deqn, init, rules, newvars},
    (* Find the equations for the initial conditions *)
     init = Select[sys, FreeQ[#, t] &];
    (* collect the differential equations *)
    deqn = Complement[sys, init];
    If[! (Equal @@ (Length /@ {deqn, init})),
      Message[NDSolve::ndnef, deqn, init];
      Return[{$Failed, $Failed, $Failed}];
      ];
    (* Find the numerical values for the initial conditions *)
    init = Solve[init, vars /. t -> t0];
    If[Not[And @@ ((VectorQ[#, NumericQ]||NumericQ[#]) & /@ (vars /. t -> t0 /. init))],
          Message[NDSolve::ndinn, init];
      Return[{$Failed, $Failed, $Failed}];
      ];
    (* Find the explicit form of the differential equations *)
    If[! FirstOrderSystemQ[deqn, vars, t],
       Message[RungeKuttaNDSolve::dord, deqn];
      Return[{$Failed, $Failed, $Failed}];
       ];
    If[! AutonomDEqnQ[deqn, vars, t],
      Message[RungeKuttaNDSolve::nauto, deqn];
      Return[{$Failed, $Failed, $Failed}];
      ];
    deqn = Solve[deqn, D[#, t] & /@ vars];
    newvars = Table[Unique[rk], {Length[vars]}];
    rules = Thread[(Head /@ vars) -> newvars];
    {deqn, init} = {deqn, init} /. rules;
    (* make pure or compiled functions from the right hand side *)
    deqn = 
      makeODEFunction[#, #[t] & /@ newvars, 
            compile] & /@ ((#'[t] &  /@ newvars) /. deqn);
    init = (#[t0] &  /@ newvars) /. init;
    {init, deqn, Reverse /@ rules}
   ]
   
(* to save time when storing the values a folded 
    list is created, so we copy always two elements *)
$interpolateOutput[t_, y_, f_, t1_, y1_, f1_, h_] :=
  ($interpolationData = {$interpolationData, {t1, y1, f1}}; True)

$contInterpolatedOutput[t_, y_, f_, t1_, y1_, f1_, h_] :=
  Block[{steps, hh, yy, ff, ii = 1},
    If[h > $RKMaxStepSize,
      steps = Floor[h/$RKMaxStepSize];
      If[Abs[hh*steps - h] <= $MachineEpsilon, steps--];
      hh = h/(steps + 1);
      $interpolationData = Nest[
            (yy = $RKInterpolationFunction[ii/(steps + 1), y, h];
              ff = $RKRighthandSide[yy];
              $FunctionCalls++;
               {#, {t + (ii++)*hh, yy, ff}}) &, $interpolationData, steps];
      ];
    $interpolationData = {$interpolationData, {t1, y1, f1}};
    True
   ]
   
$poincareOutput[t_, y_, f_, t1_, y1_, f1_, h_] :=
  Block[{s1, s2, s3, u1, u2, u3, yl = y, ym, yr = y1, i = 0},
    u1 = 0;
    u2 = 1;
    s1 = $IntersectionFunction[yl][[1]];
    s2 = $IntersectionFunction[yr][[1]];
    (* Bisection search in [t, t + h] *)
    If[Sign[s1] != Sign[s2],
      While[i++ < $PointcareIterations,
        u3 = If[s2 - s1 != 0,(* Then Regula Falsi *)
              u2 - s2*(u2 - u1)/(s2 - s1),
             (* Else Bisection *)
              (u1 + u2)/2
             ];
        If[u3 > 1 || u3 < 0, u3 = (u1 + u2)/2];
        yy = $RKInterpolationFunction[u3, y, h];
        s3 = $IntersectionFunction[yy][[1]];
        If[Sign[s1] == Sign[s3],
          s1 = s3; u1 = u3; yl = yy,
          s2 = s3; u2 = u3; yr = yy
         ];
        If[Abs[s3] <= $PoincarePrecision, Break[]];
        If[Abs[u2 - u1] <= $PoincarePrecision, Break[]];
       ];
      $interpolationData = {$interpolationData, {t + u3*h, yy}};
     ];
    True
    ]   
    
Options[RungeKuttaNDSolve] = {Compiled -> True, RKOutputFunction -> None, 
    RKStoppingFunction -> None,
    StartingStepSize -> Automatic, MaxSteps -> 10000, 
    PrecisionGoal -> 0.0000001,
    RKReturn -> RKEndPoint, Method -> Automatic, RKStatistics -> False, 
    RKInterpolationOrder -> Automatic, PSHyperplane -> 1, 
    MaxStepSize -> Infinity}

$ContinuousMethods = {rkStepVerner65, rkStep34, rkStepDormandPrince54, rkStepDormandPrince87};

(* All methods implemented *)
$RKMethods = {RKInterpolate, RKInterpolateLow, RKInterpolateMedium, 
    RKStability, RKPhaseLag}

(* We set the order when the method is selected. The $RKOrder 
    will be used when an interpolation is required. The order set here,
   will be used as option to the Interpolation[] function. *)

findRKMethod[what_, returnType_] :=
    Switch[what,
        Automatic, $RKOrder = 6; rkStepVerner65,
        RKPhaseLag, $RKOrder = 5; rkStepPhase56,
        RKStability, $RKOrder = 6; rkStepMaxStability56,
        RKInterpolateLow, $RKOrder = 3; rkStep34,
        RKInterpolateMedium, $RKOrder = 5; rkStepDormandPrince54,
        RKInterpolate, $RKOrder = 6; rkStepVerner65,
        RKInterpolateHigh, $RKOrder = 8; rkStepDormandPrince87,
    _, Message[RungeKuttaNDSolve::unmeth, $RKMethods, what]; $RKOrder = 6; 
    rkStepVerner65
    ]

findRKReturn[type_, method_] := (
    Switch[method,
      rkStepVerner65, $RKInterpolationFunction = rk65Interpolation,
      rkStep34, $RKInterpolationFunction = rk34Interpolation,
      rkStepDormandPrince54, $RKInterpolationFunction = rk54DPInterpolation,
      rkStepDormandPrince87, $RKInterpolationFunction = rk87DPInterpolation,
      _, $RKInterpolationFunction = None
      ];
    Switch[type,
        RKEndPoint, {None, Identity},
        RKInterpolation,
             $interpolationData = {};
             {If[
          MemberQ[$ContinuousMethods, 
            
            method], $contInterpolatedOutput, $interpolateOutput], \
$interpolationData &},
       RKPoincareSection,
             $interpolationData = {};
             {$poincareOutput, $interpolationData &},
       _, None]
    )

Clear[Unfold]
(* Unfold remove the folding of the interpolation output, 
  offset gives the depth of the 
    innermost element  *)
Unfold[expr_, offset_Integer] :=
  Module[{d, lst, i = 1},
       d = Depth[expr];
       lst = Table[Null, {d - offset}];
       Nest[(lst[[i++]] = Last[#]; First[#]) &, expr, d - offset];
      lst
    ]

(* modify the set {t, y, 
      y'} for the creation of the interpolation function *)
Clear[distributeTime]

distributeTime[{t_?NumericQ, v_List, dv_List}] := {t, #} & /@ 
    Transpose[{v, dv}]
distributeTime[lst : {{_?NumericQ, _List, _List} ..}] := 
  distributeTime /@ lst
    
Clear[makePSectionOut]

makePSectionOut[vars_, {_?NumericQ, val_List}] :=
    Rule @@@ Transpose[{ vars, val}]
    
Clear[makeRKReturn]

makeRKReturn[result_, RKEndPoint, replrules_, t_, t1_] :=
   Module[{resrule},  
    resrule = #[t1] & /@ (Last /@ replrules) ;
    (Rule @@@ Transpose[{resrule, #}]) & /@ result
    ]

makeRKReturn[result_, RKInterpolation, replrules_, t_, t1_] :=
  Module[{resrule, res},
       res = Reverse[Unfold[#, 3]] & /@ result;
       res = 
      Map[Interpolation[#, InterpolationOrder -> $RKOrder][t] &, 
        Transpose[distributeTime[#]] & /@ res, {2}];
       resrule = #[t] & /@ (Last /@ replrules) ;
       (Rule @@@ Transpose[{resrule, #}]) & /@ res
    ]

makeRKReturn[result_, RKPoincareSection, replrules_, _, _] :=
  Module[{res, vars},
       If[result =!= {{}},
         res = Reverse[Unfold[#, 3]] & /@ result;
         vars = Last /@ replrules;
          Map[makePSectionOut[vars, #] &, #] & /@ res,
        (* Else *)
          {}
      ]
    ]  
    
RungeKuttaNDSolve::nnstep = "StartingStepSize option must evaluate to a \
number instead of `1`."
RungeKuttaNDSolve::mxstep = "MaxSteps option must evaluate to an integer \
instead of `1`."
RungeKuttaNDSolve::prec = "PrecisionGoal option must evaluate to a floating \
point number instead of `1`."
RungeKuttaNDSolve::dord = "`1` is not a first order system of differential \
equations."
RungeKuttaNDSolve::nauto = "`1` is not a autonomous system of differential \
equations."
RungeKuttaNDSolve::unmeth = "Method must be one of `1` instead of `2`. Using \
6th order Verner method (RKInterpolate)."
RungeKuttaNDSolve::psmeth = "For Poincare sections the method must be \
RKInterpolate, or RKInterpolateLow, RKInterpolateMedium, RKInterpolateHigh instead of `1`, using RKInterpolate."

RungeKuttaNDSolve::stat = "Number of total steps: `1`, Number of fail steps: \
`2`, Number of function calls:`3`."



RungeKuttaNDSolve[eqns_, vars_List, {t_Symbol, t0_?NumericQ, t1_?NumericQ}, 
    opts___?OptionQ] :=
  Block[{$RKActiveQ = False, $RKMoreKiPresentQ = False, (* use a Block[] to free the memory of the "global values" on exit *)
         $RKMaxStepSize = Infinity, $RKRighthandSide = Null,
         $RKInterpolationFunction = None,
         $k1, $k2, $k3, $k4, $k5, $k6, $k7, $k8, $k9, $k10, $k11, $k12, $k13, $k14,
         $RKOrder = 1, $IntersectionFunction = ({1} &),
         $PoincarePrecision = 0.000001, $PointcareIterations = $IterationLimit,
         $interpolationData, $TotalSteps = 0, $FailSteps = 0, $FunctionCalls = 0,
         compile, output, h0, maxsteps, localError, init,
         funcs, replaceRules, results, returnType,
         singleRun, finishFunc = Identity, rkmethod, statistic, iporder, stopfunc},
     {compile, stopfunc, output,h0, maxsteps, localError, returnType, 
        rkmethod, statistic, $RKMaxStepSize, iporder} =
      {Compiled, RKStoppingFunction, RKOutputFunction, StartingStepSize, 
            MaxSteps, PrecisionGoal, RKReturn, Method, RKStatistics, 
            MaxStepSize, RKInterpolationOrder} /. {opts} /. Options[RungeKuttaNDSolve];
    If[returnType === RKPoincareSection,
      If[!MemberQ[{RKInterpolate, RKInterpolateLow, RKInterpolateMedium, RKInterpolateHigh, Automatic}, rkmethod],
        Message[RungeKuttaNDSolve::psmeth, rkmethod];
        rkmethod = RKInterpolation
       ];
      psfun = PSHyperplane /. {opts} /. Options[RungeKuttaNDSolve];
      $IntersectionFunction = 
      makeODEFunction[{psfun} /. a_ == b_ :> a - b, vars, compile]
     ];
    rkmethod = findRKMethod[rkmethod, returnType]; 
    If[Automatic =!= iporder && IntegerQ[iporder],
      $RKOrder = iporder
     ];
    If[! IntegerQ[maxsteps],
      Message[RungeKuttaNDSolve::mxstep, maxsteps];
      Return[$Failed];
     ];
    If[h0 === Automatic, h0 = (t1 - t0)/maxsteps];
    If[! NumericQ[h0],
      Message[RungeKuttaNDSolve::nnstep, h];
      Return[$Failed];
     ];
    If[! NumericQ[localError],
      Message[RungeKuttaNDSolve::prec, localError];
      Return [$Failed];
     ];
    $PoincarePrecision = localError;
    {output, finishFunc} = findRKReturn[returnType, rkmethod];
    {init, fun, replaceRules} = makeODESystem[eqns, vars, t, t0, compile];
    If[MemberQ[{init, fun, replaceRules}, $Failed], Return[$Failed]];
    $TotalSteps = $FailSteps = $FunctionCalls = 0;
    results = 
      Flatten[Outer[(singleRun = 
                rkDSolve[#1, #2, h0, {t0, t1}, localError, maxsteps, t, 
                  rkmethod, output, stopfunc]; finishFunc[singleRun]) &, init,
           fun, 1], 1];
    If[TrueQ[statistic],
      Message[
        RungeKuttaNDSolve::stat, $TotalSteps, $FailSteps, $FunctionCalls]
      ];
    makeRKReturn[results, returnType, replaceRules, t, t1]
    ]    

RungeKuttaNDSolve[eqns_, var_, {t_Symbol, t0_?NumericQ, t1_?NumericQ}, 
    opts___?OptionQ]:=RungeKuttaNDSolve[eqns,{var},{t,t0,t1},opts]
    
End[] (* Private *)
EndPackage[]