BeginPackage["numer`"];
Clear["numer`*"]
BCO::usage = "  BCO[a,b,c,e,order] returns order conditions of \
Numerov type methods. "
TCO::usage = "  TCO[a,b,c,e,order] returns trunc error terms of \
Numerov type methods. "
Begin["`Private`"];
Clear["numer`Private`*"];
BCO[a_, b_, c_, e_, 1] := {b . e - 1}
BCO[a_, b_, c_, e_, order_] := 
 1/S[order + 1]*
   Map[g, Map[Distribute[#] &, 
     Map[Expand[#] &, T0[a, b, c, e, order + 1], order + 1], 
     order + 1]] - 
  1/S[order + 1]*(1 + (-1)^(order + 1))/((order + 1)*order)
TCO[a_, b_, c_, e_, 1] := {b . e - 1};
TCO[a_, b_, c_, e_, order_] := 
  Module[{aa, bb, cc, ee, tr}, Off[First::normal];
   tr = Expand[
     1/Append[
        Delete[Map[First, 
          Map[Last, SSON[aa, bb, cc, ee, order]]], -1], 1] BCO[a, b, 
       c, e, order]];
   On[First::normal];
   Return[tr]];
(*--------------------------------------------------------------------------*)
RunLengthEncode[x_List] := (Through[{First, Length}[#1]] &) /@ 
   Split[x];

Combinations[list_, num_] := 
  Module[{i}, 
    Table[Map[Prepend[#, list[[i]]] &, 
      Flatten[Combinations[list, 
         num - 1][[Array[Identity, Length[list] - i + 1, i]]], 
       1], {1}], {i, 1, Length[list]}]] /; (num > 1);
Combinations[list_, 1] := Compinations[list, 1] = Map[{{#}} &, list];

Combinations2[list_, num_] := 
  Apply[Times, 
    Flatten[Combinations[list, num], 1], {1}] /; (num > 1);
Combinations2[list_, 1] := list;
(*--------------------------------------------------------------------------*)
f[x_] := 
  Apply[Times, 
    Extract[x, Map[Append[#, 1] &, Position[x, Times[_Integer, _]]]]]*
   ReplacePart[x, 1, 
    Map[Append[#, 1] &, Position[x, Times[_Integer, _]]]];
g[y_] := Map[f, y, 1];
(*--------------------------------------------------------------------------*)
T[a_, c_, e_, 0] = {e};
T[a_, c_, e_, 1] = {c};
T[a_, c_, e_, 2] = {2*a . e + c};
T[a_, c_, e_, order_] := 
  T[a, c, e, order] = 
   Module[{temp}, 
    temp = Map[Combinations2[T[a, c, e, #[[1]]], #[[2]]] &, 
      Map[RunLengthEncode[#] &, 
       IntegerPartitions[order - 2], {1}], {2}];
    temp = Map[CoverList[#] &, temp, {3}];
    temp = Apply[MyOuter, temp, {1}];
    temp = Flatten[temp, 1];
    temp = temp /. CoverList[every_] -> every;
    temp = (order - 1)*order*Map[(a . #) &, temp, {1}];
    temp = Map[(# + c) &, temp, {1}];
    temp = temp];
MyOuter[lists__] := 
  Flatten[Outer[Times, lists], Length[{lists}] - 1];
T0[a_, b_, c_, e_, 1] = {b . e - 1}; T0[a_, b_, c_, e_, 2] = {b . c};
T0[a_, b_, c_, e_, order_] := 
 T0[a, b, c, e, order] = 
  Module[{temp}, 
   temp = Map[Combinations2[T[a, c, e, #[[1]]], #[[2]]] &, 
     Map[RunLengthEncode[#] &, 
      IntegerPartitions[order - 2], {1}], {2}];
   temp = Map[CoverList[#] &, temp, {3}];
   temp = Apply[MyOuter, temp, {1}];
   temp = Flatten[temp, 1];
   temp = temp /. CoverList[every_] -> every;
   temp = Map[(b . #) &, temp, {1}]]
(*--------------------------------------------------------------------------*)
S[1] = {1}; S[2] = {1};
S[order_] := 
  S[order] = 
   Module[{temp}, 
    temp = Map[Combinations2[MapIndexed[ff, S[#[[1]]]], #[[2]]] &, 
      Map[RunLengthEncode[#] &, 
       IntegerPartitions[order - 2], {1}], {2}];
    temp = 
     temp /. {ff[a_, b_]^p_ -> Factorial[p]*a^p, ff[a_, b_] -> a};
    temp = Apply[MyOuter, temp, {1}];
    temp = Flatten[temp, 1]];
(*--------------------------------------------------------------------------*)
SSON[a_, b_, c_, e_, order_] := 
 Map[g, Map[Distribute[#] &, 
   Map[Expand[#] &, T0[a, b, c, e, order + 1], order + 1], 
   order + 1]]
End[];
EndPackage[];