(*-------------------------------------------------------------------------------------*)
NumerovLin10[fun_, x0_, xe_, ybegin0_, ybegin1_, n_] :=
Module[{b, c, a, x, h, y, m, f0, f, f1},
c = {-1, 0, -0.3662356688130395544618242275239891, 0.3662356688130395544618242275239891,
     -0.7297520822551534681767669861152164, 0.7297520822551534681767669861152164, 0, 1};
b = {0.003921568627450980392156862745098039, 0.1666666666666666666666666666666667,
     0.2409895906545774764333039199980612, 0.2409895906545774764333039199980612,
     0.08842217405130487650787255059017405, 0.08842217405130487650787255059017405,
     0.1666666666666666666666666666666667, 0.003921568627450980392156862745098039};
a = {{0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0},
     {0,-0.1160535518510525805616468787364216,0,0,0,0,0,0},
     {-0.002140044288214456548632025490273145,0.2464788082661780953924911511917838,
      0.005843352984023335056318223086056815,0,0,0,0,0},
     {0.003362375421039850515985914565199609,-0.08573899183709436900515990199364436,
      -0.01270564020140404806478750820901060,-0.003524733732252030870353899112392377,0,0,0,0},
     {13.06548607635129302282436720089670,-4,55.18969832297540758763068163968267,
      -12,-51.62403930742125773970259724921400,0,0,0},
     {-6.976608302188355093597814942916628,2.125245351726251788885634050101033,
      -29.39195588165726906650784849140723,6.613468699563125277813057280023256,
      27.62994495103068281256839451862113,-0.00009481847443571916142241442156037977,0,0},
     {1.965485508791388130945688239183845,5.563743082442459484162482974221134,
      8.905779154652039412220736923357196,-6.202589816025011714640573243915284,
     -10.00405598570782611098921466710374,0.2716380558469507983008797742568455,1/2,0}};
h = SetPrecision[(xe - x0)/n, 33];
x = x0 + Range[0, n]*h;
m = Length[ybegin0];
y = Array[Array[0 &, m] &, n + 1];
y[[1]] = ybegin0;
y[[2]] = ybegin1;
f = Array[Array[0 &, m] &, 8];
f[[2]] = fun[x[[1]], y[[1]]];
Do[f0 = f[[2]];
   f1 = fun[x[[r]], y[[r]]];
   f = Array[Array[0 &, m] &, 8];
   f[[1]] = f0;
   f[[2]] = f1;
   Do[f[[t]] = fun[x[[r]] + c[[t]]*h, (1 + c[[t]])*y[[r]]
               - c[[t]]*y[[r - 1]] + h^2*a[[t]] . f], {t, 3, 8}
     ];
   y[[r + 1]] = SetPrecision[2*y[[r]] - y[[r - 1]] + h^2*b . f, 33], {r, 2, n}
  ];
Return[{x, y}]];
(*-------------------------------------------------------------------------------------*)