%-----------------------------------------------------------------------------
NumerovSol9var[fcn_, x0_, x1_, xe_, ystart0_, ystart1_, tol_] :=
  Module[{a, b, bb, bs, c, h, x, m, y, f, f0, f1, f2, erro, k, doub,
          half, xout, yout, y0, y1, y2, stages},
   a = SetAccuracy[
   {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {1138818207635569/17888788647582928,3333378954260163/7639424778862807,
    0,0,0,0,0,0,0,0,0,0,0},
   {-738961295452201/18555377494643451,-479030696518956/4630982952248939,
    221463395534051/9040541347820963,0,0,0,0,0,0,0,0,0,0},
   {-89356055271778/8266078115467089,833921977842237/28086135942568280,
    260430744922526/8296325243696293,1917235959869338/8741740453475911,
    0,0,0,0,0,0,0,0,0},{54946321864723/7072718417900302,
    3292070461962517/15476896927194413,426310606302047/28894894451906660,
    -4872875879837389/20794769713847894,-741532042493254/7746353663357285,
    0,0,0,0,0,0,0,0},{-251224423618045/7424952508545912,
    125088230976712189/272926094612794188,218775165123335/17123639767701491,
    -666835004016636/9078203056179917,864760805293144/11070910018268767,
    960281482715860/4596707163596037,0,0,0,0,0,0,0},
   {-10576773872709440/16766773173543071,-8690145331927382/15796105007055177,
    89083185893050186/132307523842851153,12269900814562454/22837330783350247,
    -5034035877437941/6599359324070268,11389219613319343/15600320848904422,
    52072619848771/16647033931706302,0,0,0,0,0,0},{990106812662497/8608786706779221,
    47709993460294071/35784093553878749,-9736725818190157/25059812990564244,
    -160956717120047/3569163443077247,5224345334971628/20563103273326743,
    -136008468143191/15698515604364771,3475890738044597/11147604929887509,
    -1552113540014879/2714171497526320,0,0,0,0,0},
   {18416007371935979/420432758898246951,-156620231619701803/465915346979262697,
    163401025172665942/235885004367966089,12677558998648686/78953837587788515,
    -79019852266036174/470480554218988571,-38169819233982202/219967012131834185,
    -365063788309013031/906680012125429813,21925574036176765/268754978565376991,
    -41238414961974253/1762793429763497792,0,0,0,0},
   {34543668982039331/410545724822936983,286693809749635373/382376499899843436,
    -114182916006835913/337561389990306790,-25360637462765367/137962475652853258,
    -86500568247232663/235939559792002114,-19491967257348275/230984832785452168,
    11691031589066714/166076370337883255,18961131399369651/71663521282690571,
    44526216726901119/242101211697875540,-1463041441794741/310164447684055957,0,0,0},
   {124667255698865847/539737388738833342,-129716831316448379/246500686057727329,
    48349690501025998/43630540272032909,-8473022563345553/141327709915383754,
    6381224278475555/359395283617477737,-89886296068392911/654585800092602840,
    -72595035575763859/263483678412286769,38977102072170969/314686089940312541,
    -535342268702277563/1292605365671680301,-33999657291563413/201808439863022228,
    5297568021958511/740760492329053756,0,0},{28454082777317276/110034918691513293,
    680667020876240761/679623745715236211,-357139447610076715/476075397256183029,
    -111573303042065203/4510705512098184618,25086062842791575/124338736766908913,
    -7591298801595348/186554816204236381,2273910331385924/89681727067232285,
    92929059981689939/293465840017692885,21998919671668353/59906699649424436,
    -31983348834452542/232251852269372751,-40207523454342055/153739002754436812,
    -52822886315457520/176061236801621323,0}}, 33];
   bs = SetAccuracy[{-0.0017997292162861869406821038692212663,
    -0.013112505745942721318478283750136349, 0.0030595399154032869507671066168104193,
    -0.054918717796226510397919125849898161, 0.0063738383592174385411533238306056828,
    -0.040056743389166255827692962832469985, -0.0031140883015835712402402127512157241,
    0.0021309081881628802886990944047895974, 6.9715769493893019783881457848479396*10^-6,
    -0.021003288209281137177314489795767190, -0.0069913243711935013212778360310534574,
    0.0019149951354399596584310738117556711, 0.0025101438545069294825760280700159128},
    33];
   c = SetAccuracy[{-1,0,14472334024676221/23416728348467685,
     -4657513458471779/11992749807882580, 4657513458471779/11992749807882580,
     -9433332829198626/12651088817059513, 9433332829198626/12651088817059513,
     0,1,-1/2,1/2,-3/4,3/4}, 33];
   bb = SetAccuracy[{-(2663124765893345/414307129713565604),
    152575928519150547/327517848104256076, 45144011139774057/463099482327322321,
    125109351980123885/778979402784152446, 31069253439689137/277130437312338148,
    13097948064811857/111693808119155519, 5147129187992789/96921094680830912,
    0, 0, 0, 0, 0, 0}, 33];
   b = SetAccuracy[{74093748099261/21747946101134003,
       5117072112257713/13418955899996630,0, 2826332804819903/11901194670302015,
       2826332804819903/11901194670302015, 208925646744011/2633580729670742,
       208925646744011/2633580729670742,-706481309125889/32445337026320775,
       74093748099261/21747946101134003, 0, 0, 0, 0}, 33];
   h = x1 - x0; doub = 0; stages = 0;
   f0 = fcn[x0, ystart0];
   xout = Join[{x0}, {x1}, Array[0 &, 100000]];
   m = Length[ystart0];
   yout = Array[Array[0 &, m] &, 100002];
   yout[[1]] = ystart0; y0 = ystart0;
   yout[[2]] = ystart1; y1 = ystart1; k = 2;
   f = Array[Array[0 &, m] &, 13];
   f1 = fcn[x1, yout[[2]]]; x = x1; stages = 2;
   While[x < xe, f[[1]] = f0; f[[2]] = f1;
    Do[f[[o]] =
      fcn[x + c[[o]]*h, (1 + c[[o]])*y1 - c[[o]]*y0 +
        h^2*a[[o]] . f], {o, 3, 9}];
        stages = stages + 7;
    erro = 1/2*Max[Abs[h^2*(b - bb) . f]];
    If[erro < 32*tol, x = x + h; k = k + 1;               (* step accepted *)
     y2 = SetAccuracy[2*y1 - y0 + h^2*b . f, 33];
     f2 = fcn[x, y2]; stages = stages + 1;
     xout[[k]] = x; yout[[k]] = y2;
     If[erro < tol/32 && doub == 0, doub = 1; h = h + h;  (* step doubled  *)
      y1 = y2; f1 = f2, y0 = y1; f0 = f1;
      y1 = y2; f1 = f2; doub = 0],
     Do[f[[o]] =                                          (* step halved   *)
       fcn[x + c[[o]]*h, (1 + c[[o]])*y1 - c[[o]]*y0 +
                         h^2*a[[o]] . f], {o, 10, 13}];
     y0 = SetAccuracy[1/2*y1 + 1/2*y0 + h^2*bs . f, 33];
     f0 = fcn[x - h/2, y0];
     stages = stages + 5; h = h/2]];
   Return[{xout[[1 ;; k]], yout[[1 ;; k]], stages}]];

%-----------------------------------------------------------------------------