using SparseArrays
using LinearAlgebra

# DIAGONAL_SBP_D1 creates a diagonal norm SBP operator
# (D, Hinv, H, r) = diagonal_sbp_D1(p, N; xc = (-1,1))
#
# inputs:
#   p: sbp interior accuracy
#   N: finite difference grid size is N+1
#   xc: (keyword) grid span [default: (-1, 1)]
#
# outputs:
#   D:   difference operator Hinv*(-M+BS)
#   HI:  inverse of the SBP norm
#   H:   the SBP norm
#   r:   grid from xc[1] to xc[2]
#
#   References:
#   Operators for order 2, 4, 6 are from
#   @book{gustafsson2008high,
#     title={High order difference methods for time dependent PDE},
#     author={Gustafsson, Bertil},
#     year={2008},
#     publisher={Springer},
#     series={Springer Series in Computational Mathematics},
#     volume={38}
#   }
#
#   Order 8 operator is from
#   @article{MattssonSvardShoeybe2008JCP,
#     title={Stable and accurate schemes for the compressible Navier--Stokes equations},
#     author={Mattsson, Ken and Sv{\"a}rd, Magnus and Shoeybi, Mohammad},
#     journal={Journal of Computational Physics},
#     volume={227},
#     number={4},
#     pages={2293--2316},
#     year={2008},
#     doi={10.1016/j.jcp.2007.10.018}
#   }
#   with
#      x1 =  663/950
#      x2 = -229/1900
#      x3 =  415/547
#
#   Order 10 operator is from
#   @Article{MattssonAlmquist2013JCP,
#     author = {K. Mattsson and M. Almquist},
#     title = {A solution to the stability issues with block norm summation by parts operators},
#     journal = {Journal of Computational Physics},
#     volume = {253},
#     pages = {418--442},
#     year = {2013},
#     doi = {10.1016/j.jcp.2013.07.013}
#   }
#{{{
function diagonal_sbp_D1(p, N; xc = (-1, 1))

  if p == 2
    bhinv = [2]
    d  = [-1/2 0 1/2]
    bd  = [-1 1]
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49]
    d  = [1/12 -2/3 0 2/3 -1/12]
    bd = [-24/17  59/34  -4/17  -3/34  0     0;
           -1/2    0      1/2    0     0     0;
            4/43 -59/86   0     59/86 -4/43  0;
            3/98   0    -59/98   0    32/49 -4/49]
  elseif p == 6
    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    d = [-1/60 3/20 -3/4 0 3/4 -3/20 1/60];

    bd = [ -21600/13649   104009/54596   30443/81894  -33311/27298    16863/27298  -15025/163788        0           0         0;
          -104009/240260       0          -311/72078   20229/24026   -24337/48052   36661/360390        0           0         0;
           -30443/162660     311/32532       0        -11155/16266    41287/32532  -21999/54220         0           0         0;
            33311/107180  -20229/21436     485/1398        0           4147/21436    25427/321540      72/5359      0         0;
           -16863/78770    24337/31508  -41287/47262   -4147/15754        0         342523/472620   -1296/7877    144/7877    0;
            15025/525612  -36661/262806  21999/87602  -25427/262806 -342523/525612       0          32400/43801 -6480/43801 720/43801];
  elseif p == 8
    bhinv = [5080320/1498139 725760/1107307 80640/20761 725760/1304999 725760/299527 80640/103097 725760/670091 5080320/5127739];
    d = [1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 -1/280];

    bd = [    -2540160/1498139        699846290793/311403172540  -10531586157/311403172540  -145951651079/186841903524    398124597/15570158627    39152001/113858564  -80631892229/934209517620  -6230212503/311403172540        0                 0               0               0;
          -24132630717/55557028660               0                 2113176981/23016483302      5686186719/11508241651   -3408473341/138098899812  -39291999/210388330     607046586/11508241651    3460467023/483346149342        0                 0               0               0;
            3510528719/90623010660      -704392327/1294614438               0                   503511235/1294614438      354781619/2589228876       407439/3944590        -2986842/16597621        169381493/3020767022          0                 0               0               0;
          145951651079/1139279786988   -5686186719/13562854607    -1510533705/27125709214               0                6763379967/54251418428    13948923/49589962    -1603900430/40688563821   -3742312557/189879964498        0                 0               0               0;
            -398124597/21790888777      3408473341/37355809332    -1064344857/12451936444     -6763379967/12451936444             0                  763665/1198108     -1282435899/12451936444    7822226819/261490665324    -2592/299527          0               0               0;
              -1864381/23506116           13097333/58765290           -407439/19588430           -4649641/11753058          -254555/1237164               0                 5346432/9794215           -923328/9794215          3072/103097       -288/103097        0               0;
           11518841747/417855345780     -607046586/6964255763         1192698/23768791         1603900430/20892767289    1282435899/27857023052   -48117888/63658645              0                 301190400/366539777     -145152/670091      27648/670091    -2592/670091        0;
            6230212503/1065851828540   -3460467023/319755548562   -1524433437/106585182854     3742312557/106585182854  -7822226819/639511097124   58169664/487135205   -2108332800/2804873233              0               4064256/5127739  -1016064/5127739  193536/5127739  -18144/5127739];
  elseif p == 10
    bhinv = [18289152000/5261271563 1828915200/2881040311 406425600/52175551 6096384/11662993 87091200/50124587 72576000/50124587 87091200/148333439 152409600/63867949 16257024/20608675 1828915200/1704508063 18289152000/18425967263];

    d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260];

    bd = [  -1.7380923775745425e+00   2.3557601935237220e+00  -1.5328406598563976e-01  -5.7266565770416333e-01  -1.8308103515008173e-01   1.8186748267946842e-01   2.0034232582598244e-01   2.2678007363666621e-02  -1.1782459320459637e-01  -3.0591175636402144e-02   3.4890895862586133e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -4.3020203737210871e-01   0.0000000000000000e+00   1.1837297346927406e-01   3.3928601158526644e-01   1.3241927733034406e-01  -8.7495003780608913e-02  -1.1750484124279399e-01  -1.6401912273575153e-02   6.2537843443041474e-02   1.7143274696828435e-02  -1.8155585855667674e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             3.4348531361887280e-01  -1.4525207124434036e+00   0.0000000000000000e+00   2.9011513992277767e+00  -2.2419288742360557e+00  -5.4662873578741478e-01   1.2908050607446131e+00   6.1514504292452719e-02  -4.2442625460011202e-01   1.5579158905288801e-02   5.2969140277981920e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             8.6111387816878188e-02  -2.7937273515056432e-01  -1.9467880944770807e-01   0.0000000000000000e+00   2.0170150914578375e-01   2.4269917331475005e-01  -7.7261988327590472e-02   5.0649247607525059e-02  -7.4775049946661561e-03  -4.0978487203372188e-02   1.8608207238964152e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             9.1509035082611684e-02  -3.6243526359648576e-01   5.0007055839856984e-01  -6.7045605191055857e-01   0.0000000000000000e+00  -1.7807807859119628e-02   7.5000761407401195e-01  -2.2979723229714316e-01  -1.2521154324370892e-01   6.8278284106004450e-02  -4.1575927541817690e-03   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -7.5752056274147259e-02   1.9956355926115746e-01   1.0160630736447970e-01  -6.7227694623145351e-01   1.4839839882599690e-02   0.0000000000000000e+00   5.4091068834671807e-01  -1.2712520372174399e-01  -8.9292453564020990e-02   1.6181541970619609e-01  -5.4289154769785249e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -3.3838029883391296e-02   1.0867927550524317e-01  -9.7293058702223670e-02   8.6783825404790446e-02  -2.5344131542932297e-01  -2.1934035945002228e-01   0.0000000000000000e+00   2.7184438867288430e-01   1.9102691945078512e-01  -4.8646826827046824e-02  -6.2407959378425991e-03   4.6597719614658163e-04   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -1.5567948806367624e-02   6.1656604470023607e-02  -1.8844858059892756e-02  -2.3122780265804038e-01   3.1560994521078772e-01   2.0951677187991255e-01  -1.1048784865195491e+00   0.0000000000000000e+00   1.1823059621092409e+00  -5.3610400867086083e-01   1.5931375952374752e-01  -2.3673846172827626e-02   1.8939076938262100e-03   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             2.6737701764454301e-02  -7.7712278574126673e-02   4.2981266272823705e-02   1.1284579710276557e-02   5.6847566375570611e-02   4.8647834370398067e-02  -2.5665536068472994e-01  -3.9083324869946684e-01   0.0000000000000000e+00   6.5716944195909766e-01  -1.5822272208022428e-01   4.6954983762905661e-02  -7.8258306271509429e-03   6.2606645017207550e-04   0.0000000000000000e+00   0.0000000000000000e+00;
             9.4425181052687698e-03  -2.8976375375532045e-02  -2.1459742428921558e-03   8.4117843695442701e-02  -4.2165149106440383e-02  -1.1991463562335723e-01   8.8902467992349743e-02   2.4105392677971343e-01  -8.9388344421253152e-01   0.0000000000000000e+00   8.6496680152924643e-01  -2.5547312415382800e-01   6.3868281038457000e-02  -1.0644713506409501e-02   8.5157708051276015e-04   0.0000000000000000e+00;
            -9.9625965676187218e-03   2.8387641187789508e-02  -6.7495090936003027e-03  -3.5335033597892078e-02   2.3750992019053968e-03   3.7216380474824604e-02   1.0550378667904333e-02  -6.6265458456725809e-02   1.9908619649258188e-01  -8.0014409359906680e-01   0.0000000000000000e+00   8.2714572225493910e-01  -2.3632734921569687e-01   5.9081837303924217e-02  -9.8469728839873684e-03   7.8775783071898962e-04];
  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  (bm, bn) = size(bd);
  Np = N+1

  if(Np < 2*bm || Np < bn)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  @assert h > 0
  H_I = 1:Np
  H_V = ones(Np)
  H_V[1:bm] = 1 ./ bhinv[:]
  H_V[Np-bm+1:Np] = 1 ./ bhinv[end:-1:1]
  H = sparse(H_I, H_I, h * H_V)
  HI = sparse(H_I, H_I, 1 ./ (h * H_V))

  n = floor(Int64, length(d)/2);
  B_I1 = (bm+1:Np-bm) * ones(1,p+1)
  B_J1 = ones(Np-2bm,1) * (-div(p,2):div(p,2))' + B_I1
  B_V1 = ones(Np-2bm)
  B_V1 = kron(B_V1, d)

  B_I2 = (1:bm) * ones(1, bn)
  B_J2 = ones(bm) * (1:bn)'
  B_V2 = bd
  B_I3 = (Np+1) .- B_I2
  B_J3 = (Np+1) .- B_J2
  B_V3 = -B_V2
  D = sparse([B_I1[:];B_I2[:];B_I3[:]],
             [B_J1[:];B_J2[:];B_J3[:]],
             [B_V1[:];B_V2[:];B_V3[:]]/h, N+1, N+1)

  r = range(xc[1], stop=xc[2], length=N+1)

  (D, HI, H, r)
end
#}}}

#DIAGONAL_SBP_D2 creates a diagonal norm SBP operator for the 2nd drerivative
# (D, S0, SN, Hinv, H, r) = diagonal_sbp_D2(p, N; xc = (-1,1))
#
# inputs:
#   p: sbp interior accuracy
#   N: finite difference grid size is N+1
#   xc: (keyword) grid span [default: (-1, 1)]
#
# outputs:
#   D:   difference operator Hinv*(-M-S0+SN)
#   S0:  boundary derivative operator at first grid point
#   SN:  boundary derivative operator at last grid point
#   HI:  inverse of the SBP norm
#   H:   the SBP norm
#   r:   grid from xc[1] to xc[2]
#
#   References:
#   Operators for order 2, 4, 6, 8 are from
#   @book{gustafsson2008high,
#     title={High order difference methods for time dependent PDE},
#     author={Gustafsson, Bertil},
#     year={2008},
#     publisher={Springer},
#     series={Springer Series in Computational Mathematics},
#     volume={38}
#   }
#
#   Order 10 operator is from
#   @Article{MattssonAlmquist2013JCP,
#     author = {K. Mattsson and M. Almquist},
#     title = {A solution to the stability issues with block norm summation by
#     parts operators},
#     journal = {Journal of Computational Physics},
#     volume = {253},
#     pages = {418--442},
#     year = {2013},
#     doi = {10.1016/j.jcp.2013.07.013}
#   }
#{{{
function diagonal_sbp_D2(p, N; xc = (-1, 1))

  if p == 2
    bhinv = [2];

    d  = [1 -2 1];
    bd  = d;

    BS = [3/2 -2 1/2];
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49];

    d  = [-1/12 4/3 -5/2 4/3 -1/12];

    bd = [ 2    -5       4     -1       0      0;
           1    -2       1      0       0      0;
          -4/43 59/43 -110/43  59/43   -4/43   0;
          -1/49  0      59/49 -118/49  64/49  -4/49];

    BS = [11/6 -3 3/2 -1/3];

  elseif p == 6
    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    d = [1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90];

    bd = [  114170/40947   -438107/54596   336409/40947  -276997/81894     3747/13649     21035/163788      0           0         0
              6173/5860      -2066/879       3283/1758      -303/293       2111/3516       -601/4395        0           0         0
            -52391/81330    134603/32532   -21982/2711    112915/16266   -46969/16266     30409/54220       0           0         0
             68603/321540   -12423/10718   112915/32154   -75934/16077    53369/21436    -54899/160770     48/5359      0         0
             -7053/39385     86551/94524   -46969/23631    53369/15754   -87904/23631    820271/472620  -1296/7877     96/7877    0
             21035/525612   -24641/131403   30409/87602   -54899/131403  820271/525612  -117600/43801   64800/43801 -6480/43801 480/43801];

    BS = [25/12 -4 3 -4/3 1/4];

  elseif p == 8
    bhinv = [5080320/1498139 725760/1107307 80640/20761 725760/1304999 725760/299527 80640/103097 725760/670091 5080320/5127739];
    d = [-1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560]

    bd = zeros(8,12);
    bd[1, 1] =  4870382994799/1358976868290;
    bd[1, 2] =  -893640087518/75498714905  ;
    bd[1, 3] =   926594825119/60398971924  ;
    bd[1, 4] = -1315109406200/135897686829 ;
    bd[1, 5] =    39126983272/15099742981  ;
    bd[1, 6] =    12344491342/75498714905  ;
    bd[1, 7] =  -451560522577/2717953736580;

    bd[2, 1] =  333806012194/390619153855;
    bd[2, 2] = -154646272029/111605472530;
    bd[2, 3] =    1168338040/33481641759 ;
    bd[2, 4] =   82699112501/133926567036;
    bd[2, 5] =    -171562838/11160547253 ;
    bd[2, 6] =  -28244698346/167408208795;
    bd[2, 7] =   11904122576/167408208795;
    bd[2, 8] =   -2598164715/312495323084;

    bd[3, 1] =   7838984095/52731029988;
    bd[3, 2] =   1168338040/5649753213 ;
    bd[3, 3] =    -88747895/144865467  ;
    bd[3, 4] =    423587231/627750357  ;
    bd[3, 5] = -43205598281/22599012852;
    bd[3, 6] =   4876378562/1883251071 ;
    bd[3, 7] =  -5124426509/3766502142 ;
    bd[3, 8] =  10496900965/39548272491;

    bd[4, 1] =  -94978241528/828644350023;
    bd[4, 2] =   82699112501/157837019052;
    bd[4, 3] =    1270761693/13153084921 ;
    bd[4, 4] = -167389605005/118377764289;
    bd[4, 5] =   48242560214/39459254763 ;
    bd[4, 6] =  -31673996013/52612339684 ;
    bd[4, 7] =   43556319241/118377764289;
    bd[4, 8] =  -44430275135/552429566682;

    bd[5, 1] =   1455067816/21132528431;
    bd[5, 2] =   -171562838/3018932633 ;
    bd[5, 3] = -43205598281/36227191596;
    bd[5, 4] =  48242560214/9056797899 ;
    bd[5, 5] = -52276055645/6037865266 ;
    bd[5, 6] =  57521587238/9056797899 ;
    bd[5, 7] = -80321706377/36227191596;
    bd[5, 8] =   8078087158/21132528431;
    bd[5, 9] =        -1296/299527     ;

    bd[6, 1] =   10881504334/327321118845;
    bd[6, 2] =  -28244698346/140280479505;
    bd[6, 3] =    4876378562/9352031967  ;
    bd[6, 4] =  -10557998671/12469375956 ;
    bd[6, 5] =   57521587238/28056095901 ;
    bd[6, 6] = -278531401019/93520319670 ;
    bd[6, 7] =   73790130002/46760159835 ;
    bd[6, 8] = -137529995233/785570685228;
    bd[6, 9] =          2048/103097      ;
    bd[6,10] =          -144/103097      ;

    bd[7, 1] = -135555328849/8509847458140;
    bd[7, 2] =   11904122576/101307707835 ;
    bd[7, 3] =   -5124426509/13507694378  ;
    bd[7, 4] =   43556319241/60784624701  ;
    bd[7, 5] =  -80321706377/81046166268  ;
    bd[7, 6] =   73790130002/33769235945  ;
    bd[7, 7] = -950494905688/303923123505 ;
    bd[7, 8] =  239073018673/141830790969 ;
    bd[7, 9] =       -145152/670091       ;
    bd[7,10] =         18432/670091       ;
    bd[7,11] =        -1296/670091        ;


    bd[8, 1] =             0             ;
    bd[8, 2] =   -2598164715/206729925524;
    bd[8, 3] =   10496900965/155047444143;
    bd[8, 4] =  -44430275135/310094888286;
    bd[8, 5] =     425162482/2720130599  ;
    bd[8, 6] = -137529995233/620189776572;
    bd[8, 7] =  239073018673/155047444143;
    bd[8, 8] = -144648000000/51682481381 ;
    bd[8, 9] =       8128512/5127739     ;
    bd[8,10] =      -1016064/5127739     ;
    bd[8,11] =        129024/5127739     ;
    bd[8,12] =         -9072/5127739     ;

    BS = [4723/2100 -839/175 157/35 -278/105 103/140 1/175 -6/175];

  elseif p == 10
    bhinv = [18289152000/5261271563 1828915200/2881040311 406425600/52175551 6096384/11662993 87091200/50124587 72576000/50124587 87091200/148333439 152409600/63867949 16257024/20608675 1828915200/1704508063 18289152000/18425967263];

    M = zeros(11,16);

    M[ 1, 1] =  1.2056593789671863908;
    M[ 1, 2] = -1.3378814169347239658;
    M[ 1, 3] =  0.0036847309286546532061;
    M[ 1, 4] =  0.15698288365600946515;
    M[ 1, 5] = -0.0037472461482539197952;
    M[ 1, 6] = -0.0062491712449361657064;
    M[ 1, 7] = -0.029164045872729581661;
    M[ 1, 8] =  0.00054848184117832929161;
    M[ 1, 9] =  0.013613461413384884448;
    M[ 1,10] = -0.0025059220258337808220;
    M[ 1,11] = -0.00094113457993630916498;

    M[ 2, 2] =  2.1749807117105597139;
    M[ 2, 3] = -0.12369059547124894597;
    M[ 2, 4] = -0.83712574037924152603;
    M[ 2, 5] =  0.050065127254670973258;
    M[ 2, 6] =  0.0081045853127317536361;
    M[ 2, 7] =  0.097405846039248226536;
    M[ 2, 8] = -0.00068942461520402214720;
    M[ 2, 9] = -0.041326971493379188475;
    M[ 2,10] =  0.0075778529605774119402;
    M[ 2,11] =  0.0025800256160095691057;

    M[ 3, 3] =  0.18361596652499065332;
    M[ 3, 4] =  0.048289690013342693109;
    M[ 3, 5] = -0.19719621435164680412;
    M[ 3, 6] =  0.11406859029505842791;
    M[ 3, 7] = -0.029646295985488126964;
    M[ 3, 8] = -0.0016038463172861201306;
    M[ 3, 9] =  0.0032879841528337653050;
    M[ 3,10] = -0.00093242311589807387463;
    M[ 3,11] =  0.00012241332668787820533;

    M[ 4, 4] =  1.2886524606662484673;
    M[ 4, 5] = -0.14403037739488789185;
    M[ 4, 6] = -0.44846291607489015475;
    M[ 4, 7] = -0.10598334599408054277;
    M[ 4, 8] = -0.015873275740355918053;
    M[ 4, 9] =  0.073988493386459608166;
    M[ 4,10] = -0.012508848749152899785;
    M[ 4,11] = -0.0039290233894513005339;

    M[ 5, 5] =  0.51482665719685186210;
    M[ 5, 6] =  0.051199577887125103015;
    M[ 5, 7] = -0.36233561810883077365;
    M[ 5, 8] =  0.091356850268746392169;
    M[ 5, 9] =  0.0024195916108052419451;
    M[ 5,10] = -0.0018564214413731389338;
    M[ 5,11] = -0.00070192677320704413827;

    M[ 6, 6] =  0.68636003380365860083;
    M[ 6, 7] = -0.28358848290867614908;
    M[ 6, 8] = -0.13836006478253396528;
    M[ 6, 9] =  0.0076158070663111995297;
    M[ 6,10] =  0.011447010307180005164;
    M[ 6,11] = -0.0021349696610286552676;

    M[ 7, 7] =  1.5216081480839085990;
    M[ 7, 8] = -0.42653865162216293237;
    M[ 7, 9] = -0.42047484981879143123;
    M[ 7,10] =  0.019813359263872926304;
    M[ 7,11] =  0.019221397241190103344;

    M[ 8, 8] =  1.0656733504627815335;
    M[ 8, 9] = -0.66921872668484232217;
    M[ 8,10] =  0.12022033144141336599;
    M[ 8,11] = -0.030157881394591483631;

    M[ 9, 9] =  2.4064247712949611684;
    M[ 9,10] = -1.5150200315922263367;
    M[ 9,11] =  0.17373015320416595052;

    M[10,10] =  2.7682502485427255096;
    M[10,11] = -1.5975407111468405444;

    M[11,11] =  2.9033627686681129471;
    d     = [1/3150  -5/1008  5/126  -5/21  5/3 -5269/1800 5/3 -5/21 5/126 -5/1008 1/3150]
    for k = 1:5
      M[11-5+k,11 .+ (1:k)] = -d[k:-1:1];
    end

    M[1:11,1:11] = M[1:11,1:11]'+M[1:11,1:11]-diagm(0 => diag(M[1:11,1:11]));

    BS = zeros(1,16);


    BS[1:7] = -[-49/20 6 -15/2 20/3 -15/4 6/5 -1/6];
    e0 = zeros(11,1);
    e0[1] = 1;

    bd = diagm(0 => bhinv[:]) * (-M+e0*BS);

  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  (bm, bn) = size(bd);

  M = N+1;
  if(M < 2*bm || M < bn)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  @assert h > 0
  H_I = 1:M
  H_V = ones(M)
  H_V[1:bm] = bhinv[:]
  H_V[M-bm+1:M] = bhinv[end:-1:1]
  HI = sparse(H_I, H_I, H_V / h)
  H  = sparse(H_I, H_I, h ./ H_V)

  n = floor(Int64, length(d)/2);
  B_I1 = (bm+1:M-bm) * ones(1,p+1)
  B_J1 = ones(M-2bm,1) * (-div(p,2):div(p,2))' + B_I1
  B_V1 = ones(M-2bm)
  B_V1 = kron(B_V1, d)

  B_I2 = (1:bm) * ones(1, bn)
  B_J2 = ones(bm) * (1:bn)'
  B_V2 = bd
  B_I3 = (M+1) .- B_I2
  B_J3 = (M+1) .- B_J2
  B_V3 = B_V2
  D = sparse([B_I1[:];B_I2[:];B_I3[:]],
             [B_J1[:];B_J2[:];B_J3[:]],
             [B_V1[:];B_V2[:];B_V3[:]]/h^2, N+1, N+1)

  S0 = sparse(ones(length(BS)), 1:length(BS), -BS[:]/h, N+1, N+1);
  SN = sparse((N+1) * ones(length(BS)), (N+1):-1:(N+2-length(BS)),
              BS[:]/h, N+1, N+1);

  r = range(xc[1], stop=xc[2], length=N+1)

  (D, S0, SN, HI, H, r)

end
#}}}

#VARIABLE_DIAGONAL_SBP_D2 creates a diagonal norm SBP operator for the 2nd
# derivative with variable coefficients, e.g., approximates ∂/∂r( b(r) ∂u/∂r)
#
# (D, BS, Hinv, H, r) = variable_diagonal_sbp_D2(p, N, B; xc = (-1,1))
#{{{
function variable_diagonal_sbp_D2(p, N, B; xc = (-1, 1))
  r = range(xc[1], stop=xc[2], length=N+1)
  variable_diagonal_sbp_D2(p, N, B(collect(r));xc=xc)
end
function variable_diagonal_sbp_D2(p, N, B::T; xc = (-1, 1)) where T <: Number
  variable_diagonal_sbp_D2(p, N, B*ones(N+1);xc=xc)
end

function variable_diagonal_sbp_D2(p, N, B::Array{Float64,1}; xc = (-1, 1))
  @assert length(B) == N+1

  if p == 2
    bhinv = [2];

    I_M0 = [1]
    J_M0 = [1]
    V_M0 = [(B[1]+B[2])/2]

    I_MN = [N+1]
    J_MN = [N+1]
    V_MN = [(B[N]+B[N+1])/2]

    I_M = [2:N+1;
           2:N  ;
           1:N  ]
    J_M = [1:N  ;
           2:N  ;
           2:N+1]
    V_M = [-(B[1:N  ]+B[2:N+1])/2;
            (B[1:N-1]+2B[2:N]+B[3:N+1])/2;
           -(B[1:N  ]+B[2:N+1])/2];
    M = sparse([I_M0;I_M;I_MN], [J_M0;J_M;J_MN], [V_M0;V_M;V_MN])

    BS = [3/2 -2 1/2];
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49];

    BS = [11/6 -3 3/2 -1/3];

    V_M0 = zeros(6,6)

    (b1,b2,b3,b4,b5,b6,b7,b8) = B[1:8]
    V_M0[1, 1] =                (12/17)b1 + (59/192)b2 + (27010400129/345067064608)b3 + (69462376031/2070402387648)b4
    V_M0[1, 2] = V_M0[2, 1] = - (59/68)b1 - (6025413881/21126554976)b3 - (537416663/7042184992)b4
    V_M0[1, 3] = V_M0[3, 1] =   (2/17)b1 - (59/192)b2 + (213318005/16049630912)b4 + (2083938599/8024815456)b3
    V_M0[1, 4] = V_M0[4, 1] =   (3/68)b1 - (1244724001/21126554976)b3 + (752806667/21126554976)b4
    V_M0[1, 5] = V_M0[5, 1] =   (49579087/10149031312)b3 - (49579087/10149031312)b4
    V_M0[1, 6] = V_M0[6, 1] = - (1/784)b4 + (1/784)b3
    V_M0[2, 2] =                (3481/3264)b1 + (9258282831623875/7669235228057664)b3 + (236024329996203/1278205871342944)b4
    V_M0[2, 3] = V_M0[3, 2] = - (59/408)b1 - (29294615794607/29725717938208)b3 - (2944673881023/29725717938208)b4
    V_M0[2, 4] = V_M0[4, 2] = - (59/1088)b1 + (260297319232891/2556411742685888)b3 - (60834186813841/1278205871342944)b4
    V_M0[2, 5] = V_M0[5, 2] = - (1328188692663/37594290333616)b3 + (1328188692663/37594290333616)b4
    V_M0[2, 6] = V_M0[6, 2] = - (8673/2904112)b3 + (8673/2904112)b4
    V_M0[3, 3] =                (1/51)b1 + (59/192)b2 + (13777050223300597/26218083221499456)b4 + (564461/13384296)b5 + (378288882302546512209/270764341349677687456)b3
    V_M0[3, 4] = V_M0[4, 3] =   (1/136)b1 - (125059/743572)b5 - (4836340090442187227/5525802884687299744)b3 - (17220493277981/89177153814624)b4
    V_M0[3, 5] = V_M0[5, 3] = - (10532412077335/42840005263888)b4 + (1613976761032884305/7963657098519931984)b3 + (564461/4461432)b5
    V_M0[3, 6] = V_M0[6, 3] = - (960119/1280713392)b4 - (3391/6692148)b5 + (33235054191/26452850508784)b3
    V_M0[4, 4] =                (3/1088)b1 + (507284006600757858213/475219048083107777984)b3 + (1869103/2230716)b5 + (1/24)b6 + (1950062198436997/3834617614028832)b4
    V_M0[4, 5] = V_M0[5, 4] = - (4959271814984644613/20965546238960637264)b3 - (1/6)b6 - (15998714909649/37594290333616)b4 - (375177/743572)b5
    V_M0[4, 6] = V_M0[6, 4] = - (368395/2230716)b5 + (752806667/539854092016)b3 + (1063649/8712336)b4 + (1/8)b6
    V_M0[5, 5] =                (8386761355510099813/128413970713633903242)b3 + (2224717261773437/2763180339520776)b4 + (5/6)b6 + (1/24)b7 + (280535/371786)b5
    V_M0[5, 6] = V_M0[6, 5] = - (35039615/213452232)b4 - (1/6)b7 - (13091810925/13226425254392)b3 - (1118749/2230716)b5 - (1/2)b6
    V_M0[6, 6] =                (3290636/80044587)b4 + (5580181/6692148)b5 + (5/6)b7 + (1/24)b8 + (660204843/13226425254392)b3 + (3/4)b6
    I_M0 = (1:6) * ones(1,6)
    J_M0 = I_M0'

    V_MN = zeros(6,6)
    (b1,b2,b3,b4,b5,b6,b7,b8) = B[N+1:-1:N-6]
    V_MN[6, 6] =                (12/17)b1 + (59/192)b2 + (27010400129/345067064608)b3 + (69462376031/2070402387648)b4
    V_MN[6, 5] = V_MN[5, 6] = - (59/68)b1 - (6025413881/21126554976)b3 - (537416663/7042184992)b4
    V_MN[6, 4] = V_MN[4, 6] =   (2/17)b1 - (59/192)b2 + (213318005/16049630912)b4 + (2083938599/8024815456)b3
    V_MN[6, 3] = V_MN[3, 6] =   (3/68)b1 - (1244724001/21126554976)b3 + (752806667/21126554976)b4
    V_MN[6, 2] = V_MN[2, 6] =   (49579087/10149031312)b3 - (49579087/10149031312)b4
    V_MN[6, 1] = V_MN[1, 6] = - (1/784)b4 + (1/784)b3
    V_MN[5, 5] =                (3481/3264)b1 + (9258282831623875/7669235228057664)b3 + (236024329996203/1278205871342944)b4
    V_MN[5, 4] = V_MN[4, 5] = - (59/408)b1 - (29294615794607/29725717938208)b3 - (2944673881023/29725717938208)b4
    V_MN[5, 3] = V_MN[3, 5] = - (59/1088)b1 + (260297319232891/2556411742685888)b3 - (60834186813841/1278205871342944)b4
    V_MN[5, 2] = V_MN[2, 5] = - (1328188692663/37594290333616)b3 + (1328188692663/37594290333616)b4
    V_MN[5, 1] = V_MN[1, 5] = - (8673/2904112)b3 + (8673/2904112)b4
    V_MN[4, 4] =                (1/51)b1 + (59/192)b2 + (13777050223300597/26218083221499456)b4 + (564461/13384296)b5 + (378288882302546512209/270764341349677687456)b3
    V_MN[4, 3] = V_MN[3, 4] =   (1/136)b1 - (125059/743572)b5 - (4836340090442187227/5525802884687299744)b3 - (17220493277981/89177153814624)b4
    V_MN[4, 2] = V_MN[2, 4] = - (10532412077335/42840005263888)b4 + (1613976761032884305/7963657098519931984)b3 + (564461/4461432)b5
    V_MN[4, 1] = V_MN[1, 4] = - (960119/1280713392)b4 - (3391/6692148)b5 + (33235054191/26452850508784)b3
    V_MN[3, 3] =                (3/1088)b1 + (507284006600757858213/475219048083107777984)b3 + (1869103/2230716)b5 + (1/24)b6 + (1950062198436997/3834617614028832)b4
    V_MN[3, 2] = V_MN[2, 3] = - (4959271814984644613/20965546238960637264)b3 - (1/6)b6 - (15998714909649/37594290333616)b4 - (375177/743572)b5
    V_MN[3, 1] = V_MN[1, 3] = - (368395/2230716)b5 + (752806667/539854092016)b3 + (1063649/8712336)b4 + (1/8)b6
    V_MN[2, 2] =                (8386761355510099813/128413970713633903242)b3 + (2224717261773437/2763180339520776)b4 + (5/6)b6 + (1/24)b7 + (280535/371786)b5
    V_MN[2, 1] = V_MN[1, 2] = - (35039615/213452232)b4 - (1/6)b7 - (13091810925/13226425254392)b3 - (1118749/2230716)b5 - (1/2)b6
    V_MN[1, 1] =                (3290636/80044587)b4 + (5580181/6692148)b5 + (5/6)b7 + (1/24)b8 + (660204843/13226425254392)b3 + (3/4)b6

    I_MN = (N-4:N+1) * ones(1,6)
    J_MN = I_MN'

    I_M = [7:N-3;
           7:N-4;
           7:N-5;
           6:N-5;
           5:N-5]
    J_M = [5:N-5;
           6:N-5;
           7:N-5;
           7:N-4;
           7:N-3]
    V_M = [(+(1/8)B[(5:N-5).+2] - (1/6)B[(5:N-5).+1] + (1/8)B[(5:N-5)]);
           (-(1/6)B[(6:N-5).+2] - (1/2)B[(6:N-5).+1] - (1/2)B[(6:N-5)] - (1/6)B[(6:N-5).-1]);
           (+(1/24)B[(7:N-5).+2] + (5/6)B[(7:N-5).+1] + (3/4)B[(7:N-5)] + (5/6)B[(7:N-5).-1] + (1/24)B[(7:N-5).-2]);
           (-(1/6)B[(7:N-4).+1] - (1/2)B[(7:N-4)] - (1/2)B[(7:N-4).-1] - (1/6)B[(7:N-4).-2]);
           (+(1/8)B[(7:N-3)] -(1/6)B[(7:N-3).-1] +(1/8)B[(7:N-3).-2]);
          ]

    M = sparse([I_M0[:];I_M;I_MN[:]], [J_M0[:];J_M;J_MN[:]],
               [V_M0[:];V_M;V_MN[:]])

  elseif p == 6

    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    BS = [25/12 -4 3 -4/3 1/4];

    V_M0 = zeros(9,9)

    (b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12) = B[1:12]
    V_M0[1, 1] = 0.79126675946955820939b1 + 0.29684720906380007429b2 + 0.0031855190887964290152b3 + 0.016324040425909519534b4 + 0.031603022440944150877b5 + 0.031679647480161052996b6 + 0.031485777339472539205b7
    V_M0[1, 2] = V_M0[2, 1] = -1.0166893393503381444b1 - 0.028456273704916113690b3 - 0.041280298383492988198b4 - 0.13922814516201405075b5 -0.11957773256112017666b6 -0.11942677565293334109b7
    V_M0[1, 3] = V_M0[3, 1] = 0.070756429372437150463b1 - 0.18454761060241510503b2 - 0.043641631471118923470b4 + 0.24323679072077324609b5 + 0.15821270735372154440b6 + 0.16023485783647863076b7
    V_M0[1, 4] = V_M0[4, 1] = 0.22519915328913532127b1 - 0.16627487110970548953b2 + 0.027105309616486712977b3 - 0.19166461859684399091b5 - 0.076841171601990145944b6 - 0.082195869498316975759b7
    V_M0[1, 5] = V_M0[5, 1] = -0.052244034642020563167b1 + 0.044400639485098762210b2 - 0.0010239765473093878745b3 + 0.074034846453161740905b4 + 0.012416255689984968954b6  +0.071886528478926012827b5 + 0.013793629971047355034b7
    V_M0[1, 6] = V_M0[6, 1] = -0.018288968138771973527b1 +0.0095746331632217580607b2 - 0.00081057845305764042779b3 - 0.0073488455877755196984b4 + 0.010636019497239069970b5 - 0.013159670383826183824b6 - 0.021179364788387535246b7
    V_M0[1, 7] = V_M0[7, 1] = 0.0019118885633161709274b4 - 0.040681303555291499361b5 + 0.013196749810737491670b6 + 0.025572665181237836763b7
    V_M0[1, 8] = V_M0[8, 1] = 0.015596528711367857640b5 -0.0064861841573315378995b6 -0.0091103445540363197401b7
    V_M0[1, 9] = V_M0[9, 1] = 0.00055939836966298630593b6 - 0.0013848225351007963723b5 + 0.00082542416543781006633b7
    V_M0[2, 2] =              1.3063321571116676286b1 +0.25420017604573457435b3 + 0.10438978280925626095b4 + 0.66723280210321129509b5 + 0.46818193597227494411b6 + 0.46764154101958369201b7
    V_M0[2, 3] = V_M0[3, 2] = -0.090914102699924646049b1 +0.11036113131714764253b4 -1.2903975449975188870b5 -0.66396052487350447871b6 -0.66159744640052061842b7
    V_M0[2, 4] = V_M0[4, 2] = -0.28935573956534316666b1 -0.24213200040645927216b3 +1.1876702550280310277b5 +0.39565981499041363328b6 +0.38600489217558000007b7
    V_M0[2, 5] = V_M0[5, 2] = 0.067127744758037639890b1 +0.0091471926820756301800b3 -0.18721961430038080217b4 -0.13193585588531745301b6 -0.48715757368119118874b5 -0.10475163122754481381b7
    V_M0[2, 6] = V_M0[6, 2] = 0.023499279745900688694b1 +0.0072409053835651813164b3 +0.018583789963916794487b4-0.092896161339386761743b5 +0.12235132704188076670b6 +0.11135203204362950339b7
    V_M0[2, 7] = V_M0[7, 2] = -0.0048347914064469075906b4 +0.23106838326878204031b5 -0.10807741421960079917b6-0.11815617764273433354b7
    V_M0[2, 8] = V_M0[8, 2] = -0.083681414344034553537b5 +0.040934994667670546616b6 +0.042746419676364006921b7
    V_M0[2, 9] = V_M0[9, 2] = -0.0035765451326969831434b6 +0.0073893991241210786821b5 -0.0038128539914240955387b7
    V_M0[3, 3] =              0.0063271611471368738078b1 +0.11473182007158685275b2 +0.11667405542796800075b4+2.7666108082854440372b5 +1.0709206899608171042b6 +1.0131613910329730572b7
    V_M0[3, 4] = V_M0[4, 3] = 0.020137694138847972466b1 +0.10337179946308864017b2 -2.9132216211517427243b5-0.87558073434822622598b6 -0.69099571834888124265b7
    V_M0[3, 5] = V_M0[5, 3] = -0.0046717510915754628683b1 -0.027603533656377128278b2 -0.19792902986208699745b4 +0.54029853383734330523b6 +1.2391775930319110779b5 +0.26280380502473582273b7
    V_M0[3, 6] = V_M0[6, 3] = -0.0016354308669218878195b1 -0.0059524752758832596197b2 +0.019646827777442752194b4+0.32366400126390466006b5 -0.46595166932288709739b6 -0.22172727209417368594b7
    V_M0[3, 7] = V_M0[7, 3] = -0.0051113531893524745496b4 -0.53558781637747543460b5 +0.33283351044897389336b6+0.20786565911785401579b7
    V_M0[3, 8] = V_M0[8, 3] = 0.18243281741342895622b5 -0.10598160301968184459b6 -0.076451214393747111630b7
    V_M0[3, 9] = V_M0[9, 3] = 0.0092090899634437994856b6 -0.015915028188724931671b5 +0.0067059382252811321853b7
    V_M0[4, 4] =              0.064092997759871869867b1 +0.093136576388046999489b2 +0.23063676246347492291b3+3.6894403082837166203b5 +1.1905503386876088738b6 +0.59124795468888565194b7
    V_M0[4, 5] = V_M0[5, 4] = -0.014868958192656041286b1 -0.024870405993901607642b2 -0.0087129289077117541871b3 -1.2635078373718242057b6 -0.30583173978439973269b7 -1.4706919260458029548b5
    V_M0[4, 6] = V_M0[6, 4] = -0.0052051474298559556576b1 -0.0053630987475285424890b2 -0.0068971427657906095463b3-0.78575245216674501017b5 +0.22911480054237346001b7 +0.99770643562927505292b6
    V_M0[4, 7] = V_M0[7, 4] = 0.66972974880676622652b5 -0.50132473560721279390b6 -0.17951612431066454373b7
    V_M0[4, 8] = V_M0[8, 4] = -0.20229090601117515652b5 +0.14534218580636584986b6 +0.056948720204809306656b7
    V_M0[4, 9] = V_M0[9, 4] = -0.012004296184410038337b6 -0.0047769156693859238415b7 +0.016781211853795962179b5
    V_M0[5, 5] =              0.0034494550959102336252b1 +0.0066411834994278261016b2 +0.00032915450832718628585b3+0.33577217075764772000b4 +2.0964133295790264390b6 +0.23173232041831268550b7 +0.0061078257643682645765b8 +0.71091258506833766956b5
    V_M0[5, 6] = V_M0[6, 5] = 0.0012075440723041938061b1 +0.0014321166657521476075b2 +0.00026055826461832559573b3-0.033329411132516353908b4 -0.28082416973855326836b7 -0.027209080835250836084b8+0.10458654356829219874b5 -1.3484369866671155432b6
    V_M0[5, 7] = V_M0[7, 5] = 0.0086710380841746926251b4 +0.17360734113554285637b6 +0.053313621252876254126b8 -0.24249352624045263018b5 +0.15690152576785882706b7
    V_M0[5, 8] = V_M0[8, 5] = -0.086316839802171222760b6 +0.026988423604709992435b7 +0.080981941477156510853b5 -0.032764636390806391639b8
    V_M0[5, 9] = V_M0[9, 5] = 0.0074620594845308550733b6 -0.00081216403616686789496b7 +0.00055227020881270902093b8-0.0072021656571766961993b5
    V_M0[6, 6] =              0.00042272261734493450425b1 +0.00030882419443789644048b2 +0.00020625757066474306202b3+0.0033083434042009682567b4 +0.58280470164050018158b5 +0.80541742203662154736b7+0.13383632334100334433b8 +0.0055555555555555555556b9 +1.1903620718618930511b6
    V_M0[6, 7] = V_M0[7, 6] = -0.00086070442526864133026b4 -0.17480747086739049893b5 -0.31325448501150501650b8 -0.025000000000000000000b9 -0.31691663053104292713b7 -0.66916070916479291611b6
    V_M0[6, 8] = V_M0[8, 6] = 0.033546617916933521087b5 -0.33436200223869714050b7 +0.050000000000000000000b9+0.21697906098076027508b6 +0.18383632334100334433b8
    V_M0[6, 9] = V_M0[9, 6] = 0.029125184768230046430b7 +0.022790919164749163916b8 -0.030689859975187405305b6 -0.0017817995133473605962b5 -0.030555555555555555556b9
    V_M0[7, 7] =              0.00022392237357715991790b4 +0.12754377854309566738b5 +1.0116994839296081646b6 +0.96988172751725752475b8 +0.12500000000000000000b9 +0.0055555555555555555556b10 +0.48231775430312815001b7
    V_M0[7, 8] = V_M0[8, 7] = -0.037841139730330129499b5 -0.29975568851348273616b6 -0.30000000000000000000b9 - 0.025000000000000000000b10 -0.39914868674468211784b7 -0.43825448501150501650b8
    V_M0[7, 9] = V_M0[9, 7] = 0.046981462180226839339b6 -0.29668637874712374587b8 +0.050000000000000000000b10 +0.17163557041460064817b7 +0.0030693461522962583624b5 +0.17500000000000000000b9
    V_M0[8, 8] =              0.012303289427168044554b5 +0.11836475296458983325b6 +0.94105118982279433342b7 +0.95000000000000000000b9 +0.12500000000000000000b10 +0.0055555555555555555556b11 +0.56994743445211445545b8
    V_M0[8, 9] = V_M0[9, 8] = -0.023080678926719163396b6 -0.29866250537751494972b7 -0.30000000000000000000b10 -0.025000000000000000000b11 -0.0010477348605150508026b5 -0.42720908083525083608b8 - 0.42500000000000000000b9
    V_M0[9, 9] =              0.0051393702211491099770b6 +0.12477232150094220014b7 +0.95055227020881270902b8 +0.95000000000000000000b10 +0.12500000000000000000b11 +0.0055555555555555555556b12 +0.000091593624651536418269b5 +0.56111111111111111111b9
    I_M0 = (1:9) * ones(1,9)
    J_M0 = I_M0'

    V_MN = zeros(9,9)
    (b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12) = B[N+1:-1:N-10]
    V_MN[9, 9] =              0.79126675946955820939b1 + 0.29684720906380007429b2 + 0.0031855190887964290152b3 + 0.016324040425909519534b4 + 0.031603022440944150877b5 + 0.031679647480161052996b6 + 0.031485777339472539205b7
    V_MN[9, 8] = V_MN[8, 9] = -1.0166893393503381444b1 - 0.028456273704916113690b3 - 0.041280298383492988198b4 - 0.13922814516201405075b5 -0.11957773256112017666b6 -0.11942677565293334109b7
    V_MN[9, 7] = V_MN[7, 9] = 0.070756429372437150463b1 - 0.18454761060241510503b2 - 0.043641631471118923470b4 + 0.24323679072077324609b5 + 0.15821270735372154440b6 + 0.16023485783647863076b7
    V_MN[9, 6] = V_MN[6, 9] = 0.22519915328913532127b1 - 0.16627487110970548953b2 + 0.027105309616486712977b3 - 0.19166461859684399091b5 - 0.076841171601990145944b6 - 0.082195869498316975759b7
    V_MN[9, 5] = V_MN[5, 9] = -0.052244034642020563167b1 + 0.044400639485098762210b2 - 0.0010239765473093878745b3 + 0.074034846453161740905b4 + 0.012416255689984968954b6  +0.071886528478926012827b5 + 0.013793629971047355034b7
    V_MN[9, 4] = V_MN[4, 9] = -0.018288968138771973527b1 +0.0095746331632217580607b2 -0.00081057845305764042779b3 -0.0073488455877755196984b4 +0.010636019497239069970b5 -0.013159670383826183824b6-0.021179364788387535246b7
    V_MN[9, 3] = V_MN[3, 9] = 0.0019118885633161709274b4 -0.040681303555291499361b5 +0.013196749810737491670b6 +0.025572665181237836763b7
    V_MN[9, 2] = V_MN[2, 9] = 0.015596528711367857640b5 -0.0064861841573315378995b6 -0.0091103445540363197401b7
    V_MN[9, 1] = V_MN[1, 9] = 0.00055939836966298630593b6 -0.0013848225351007963723b5 +0.00082542416543781006633b7
    V_MN[8, 8] =              1.3063321571116676286b1 +0.25420017604573457435b3 +0.10438978280925626095b4 +0.66723280210321129509b5 +0.46818193597227494411b6 +0.46764154101958369201b7
    V_MN[8, 7] = V_MN[7, 8] = -0.090914102699924646049b1 +0.11036113131714764253b4 -1.2903975449975188870b5 -0.66396052487350447871b6 -0.66159744640052061842b7
    V_MN[8, 6] = V_MN[6, 8] = -0.28935573956534316666b1 -0.24213200040645927216b3 +1.1876702550280310277b5 +0.39565981499041363328b6 +0.38600489217558000007b7
    V_MN[8, 5] = V_MN[5, 8] = 0.067127744758037639890b1 +0.0091471926820756301800b3 -0.18721961430038080217b4 -0.13193585588531745301b6 -0.48715757368119118874b5 -0.10475163122754481381b7
    V_MN[8, 4] = V_MN[4, 8] = 0.023499279745900688694b1 +0.0072409053835651813164b3 +0.018583789963916794487b4-0.092896161339386761743b5 +0.12235132704188076670b6 +0.11135203204362950339b7
    V_MN[8, 3] = V_MN[3, 8] = -0.0048347914064469075906b4 +0.23106838326878204031b5 -0.10807741421960079917b6-0.11815617764273433354b7
    V_MN[8, 2] = V_MN[2, 8] = -0.083681414344034553537b5 +0.040934994667670546616b6 +0.042746419676364006921b7
    V_MN[8, 1] = V_MN[1, 8] = -0.0035765451326969831434b6 +0.0073893991241210786821b5 -0.0038128539914240955387b7
    V_MN[7, 7] =               0.0063271611471368738078b1 +0.11473182007158685275b2 +0.11667405542796800075b4+2.7666108082854440372b5 +1.0709206899608171042b6 +1.0131613910329730572b7
    V_MN[7, 6] = V_MN[6, 7] = 0.020137694138847972466b1 +0.10337179946308864017b2 -2.9132216211517427243b5-0.87558073434822622598b6 -0.69099571834888124265b7
    V_MN[7, 5] = V_MN[5, 7] = -0.0046717510915754628683b1 -0.027603533656377128278b2 -0.19792902986208699745b4 +0.54029853383734330523b6 +1.2391775930319110779b5 +0.26280380502473582273b7
    V_MN[7, 4] = V_MN[4, 7] = -0.0016354308669218878195b1 -0.0059524752758832596197b2 +0.019646827777442752194b4+0.32366400126390466006b5 -0.46595166932288709739b6 -0.22172727209417368594b7
    V_MN[7, 3] = V_MN[3, 7] = -0.0051113531893524745496b4 -0.53558781637747543460b5 +0.33283351044897389336b6+0.20786565911785401579b7
    V_MN[7, 2] = V_MN[2, 7] = 0.18243281741342895622b5 -0.10598160301968184459b6 -0.076451214393747111630b7
    V_MN[7, 1] = V_MN[1, 7] = 0.0092090899634437994856b6 -0.015915028188724931671b5 +0.0067059382252811321853b7
    V_MN[6, 6] =              0.064092997759871869867b1 +0.093136576388046999489b2 +0.23063676246347492291b3+3.6894403082837166203b5 +1.1905503386876088738b6 +0.59124795468888565194b7
    V_MN[6, 5] = V_MN[5, 6] = -0.014868958192656041286b1 -0.024870405993901607642b2 -0.0087129289077117541871b3 -1.2635078373718242057b6 -0.30583173978439973269b7 -1.4706919260458029548b5
    V_MN[6, 4] = V_MN[4, 6] = -0.0052051474298559556576b1 -0.0053630987475285424890b2 -0.0068971427657906095463b3-0.78575245216674501017b5 +0.22911480054237346001b7 +0.99770643562927505292b6
    V_MN[6, 3] = V_MN[3, 6] = 0.66972974880676622652b5 -0.50132473560721279390b6 -0.17951612431066454373b7
    V_MN[6, 2] = V_MN[2, 6] = -0.20229090601117515652b5 +0.14534218580636584986b6 +0.056948720204809306656b7
    V_MN[6, 1] = V_MN[1, 6] = -0.012004296184410038337b6 -0.0047769156693859238415b7 +0.016781211853795962179b5
    V_MN[5, 5] =              0.0034494550959102336252b1 +0.0066411834994278261016b2 +0.00032915450832718628585b3+0.33577217075764772000b4 +2.0964133295790264390b6 +0.23173232041831268550b7 +0.0061078257643682645765b8 +0.71091258506833766956b5
    V_MN[5, 4] = V_MN[4, 5] = 0.0012075440723041938061b1 +0.0014321166657521476075b2 +0.00026055826461832559573b3-0.033329411132516353908b4 -0.28082416973855326836b7 -0.027209080835250836084b8+0.10458654356829219874b5 -1.3484369866671155432b6
    V_MN[5, 3] = V_MN[3, 5] = 0.0086710380841746926251b4 +0.17360734113554285637b6 +0.053313621252876254126b8 -0.24249352624045263018b5 +0.15690152576785882706b7
    V_MN[5, 2] = V_MN[2, 5] = -0.086316839802171222760b6 +0.026988423604709992435b7 +0.080981941477156510853b5 -0.032764636390806391639b8
    V_MN[5, 1] = V_MN[1, 5] = 0.0074620594845308550733b6 -0.00081216403616686789496b7 +0.00055227020881270902093b8-0.0072021656571766961993b5
    V_MN[4, 4] =              0.00042272261734493450425b1 +0.00030882419443789644048b2 +0.00020625757066474306202b3+0.0033083434042009682567b4 +0.58280470164050018158b5 +0.80541742203662154736b7+0.13383632334100334433b8 +0.0055555555555555555556b9 +1.1903620718618930511b6
    V_MN[4, 3] = V_MN[3, 4] = -0.00086070442526864133026b4 -0.17480747086739049893b5 -0.31325448501150501650b8 -0.025000000000000000000b9 -0.31691663053104292713b7 -0.66916070916479291611b6
    V_MN[4, 2] = V_MN[2, 4] = 0.033546617916933521087b5 -0.33436200223869714050b7 +0.050000000000000000000b9+0.21697906098076027508b6 +0.18383632334100334433b8
    V_MN[4, 1] = V_MN[1, 4] = 0.029125184768230046430b7 +0.022790919164749163916b8 -0.030689859975187405305b6 -0.0017817995133473605962b5 -0.030555555555555555556b9
    V_MN[3, 3] =              0.00022392237357715991790b4 +0.12754377854309566738b5 +1.0116994839296081646b6 +0.96988172751725752475b8 +0.12500000000000000000b9 +0.0055555555555555555556b10 +0.48231775430312815001b7
    V_MN[3, 2] = V_MN[2, 3] = -0.037841139730330129499b5 -0.29975568851348273616b6 -0.30000000000000000000b9 - 0.025000000000000000000b10 -0.39914868674468211784b7 -0.43825448501150501650b8
    V_MN[3, 1] = V_MN[1, 3] = 0.046981462180226839339b6 -0.29668637874712374587b8 +0.050000000000000000000b10 +0.17163557041460064817b7 +0.0030693461522962583624b5 +0.17500000000000000000b9
    V_MN[2, 2] =               0.012303289427168044554b5 +0.11836475296458983325b6 +0.94105118982279433342b7 +0.95000000000000000000b9 +0.12500000000000000000b10 +0.0055555555555555555556b11 +0.56994743445211445545b8
    V_MN[2, 1] = V_MN[1, 2] = -0.023080678926719163396b6 -0.29866250537751494972b7 -0.30000000000000000000b10 -0.025000000000000000000b11 -0.0010477348605150508026b5 -0.42720908083525083608b8 - 0.42500000000000000000b9
    V_MN[1, 1] =              0.0051393702211491099770b6 +0.12477232150094220014b7 +0.95055227020881270902b8 +0.95000000000000000000b10 +0.12500000000000000000b11 +0.0055555555555555555556b12 +0.000091593624651536418269b5 +0.56111111111111111111b9

    I_MN = (N-7:N+1) * ones(1,9)
    J_MN = I_MN'

    I_M = [10:N-5;
        10:N-6;
	10:N-7;
	10:N-8;
	9:N-8;
	8:N-8;
	7:N-8]
    J_M = [7:N-8;
        8:N-8;
	9:N-8;
	10:N-8;
	10:N-7;
	10:N-6;
	10:N-5]
    V_M =  [(-(11/360)B[(10:N-5).-3] + (1/40)B[(10:N-5).-2] + (1/40)B[(10:N-5).-1] - (11/360)B[(10:N-5)]);
    	    (+(1/20)B[(10:N-6).-3] + (7/40)B[(10:N-6).-2] - (3/10)B[(10:N-6).-1] + (7/40)B[(10:N-6)] + (1/20)B[(10:N-6).+1]);
            (-(1/40)B[(10:N-7).-3] -(3/10)B[(10:N-7).-2] - (17/40)B[(10:N-7).-1] - (17/40)B[(10:N-7)] - (3/10)B[(10:N-7).+1] - (1/40)B[(10:N-7).+2]);
	    (+(1/180)B[(10:N-8).-3] +(1/8)B[(10:N-8).-2] + (19/20)B[(10:N-8).-1] + (101/180)B[(10:N-8)] + (19/20)B[(10:N-8).+1] + (1/8)B[(10:N-8).+2]+ (1/180)B[(10:N-8).+3]);
	    (-(1/40)B[(9:N-8).-2] - (3/10)B[(9:N-8).-1] - (17/40)B[(9:N-8)] - (17/40)B[(9:N-8).+1] - (3/10)B[(9:N-8).+2] - (1/40)B[(9:N-8).+3]);
	    (+(1/20)B[(8:N-8).-1] + (7/40)B[(8:N-8)] - (3/10)B[(8:N-8).+1] + (7/40)B[(8:N-8).+2]+ (1/20)B[(8:N-8).+3] );
	    (-(11/360)B[(7:N-8)] + (1/40)B[(7:N-8).+1] + (1/40)B[(7:N-8).+2] - (11/360)B[(7:N-8).+3]);
          ]

    M = sparse([I_M0[:];I_M[:];I_MN[:]],
               [J_M0[:];J_M[:];J_MN[:]],
               [V_M0[:];V_M[:];V_MN[:]])


  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  bm = length(bhinv)
  Np = N+1

  if(Np < 2*bm)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  M = M/h
  @assert h > 0
  H_I = 1:Np
  H_V = ones(Np)
  H_V[1:bm] = bhinv[:]
  H_V[Np-bm+1:Np] = bhinv[end:-1:1]
  HI = sparse(H_I, H_I, H_V / h)
  H  = sparse(H_I, H_I, h ./ H_V)

  S0 = sparse(ones(length(BS)), 1:length(BS), -B[1]*BS[:]/h, N+1, N+1);
  SN = sparse((N+1) * ones(length(BS)), (N+1):-1:(N+2-length(BS)),
              B[N+1]*BS[:]/h, N+1, N+1);
  D = HI * (-M + SN - S0)

  r = range(xc[1], stop=xc[2], length=N+1)

  BS = SN - S0
  (D, BS, HI, H, r)
end
#}}}

nothing
