--- a/+sbp/+implementations/d4_compatible_6.m	Fri Sep 09 11:03:13 2016 +0200
+++ b/+sbp/+implementations/d4_compatible_6.m	Fri Sep 09 14:53:41 2016 +0200
@@ -31,11 +31,11 @@
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
     end
 
-    H=diag(ones(m,1),0);
+    H=speye(m,m);
     H_U=[0.7493827e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5534051e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.104561e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.260503e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.43237e5 / 0.103680e6 0 0 0; 0 0 0 0 0 0.514081e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3356179e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25631027e8 / 0.25401600e8;];
 
     H(1:8,1:8)=H_U;
-    H(m-7:m,m-7:m)=fliplr(flipud(H_U));
+    H(m-7:m,m-7:m)=rot90(H_U,2);
     H=H*h;
     HI=inv(H);
 
@@ -43,18 +43,22 @@
     % First derivative SBP operator, 3rd order accurate at first 8 boundary points
 
     q3=1/60;q2=-3/20;q1=3/4;
-    Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+%     Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
+    stencil = [-q3,-q2,-q1,0,q1,q2,q3];
+    d = (length(stencil)-1)/2;
+    diags = -d:d;
+    Q = stripeMatrix(stencil, diags, m);
 
     Q_U = [0 0.26431903e8 / 0.43545600e8 0.39791489e8 / 0.152409600e9 -0.12747751e8 / 0.16934400e8 0.76099447e8 / 0.152409600e9 -0.1397443e7 / 0.12192768e8 0 0; -0.26431903e8 / 0.43545600e8 0 -0.13847476213e11 / 0.19559232000e11 0.35844843977e11 / 0.11735539200e11 -0.63413503537e11 / 0.23471078400e11 0.4764412871e10 / 0.3911846400e10 -0.1668252557e10 / 0.5867769600e10 0.842644697e9 / 0.29338848000e11; -0.39791489e8 / 0.152409600e9 0.13847476213e11 / 0.19559232000e11 0 -0.73834802771e11 / 0.23471078400e11 0.1802732209e10 / 0.325987200e9 -0.65514173e8 / 0.16299360e8 0.79341409141e11 / 0.58677696000e11 -0.1282384321e10 / 0.7823692800e10; 0.12747751e8 / 0.16934400e8 -0.35844843977e11 / 0.11735539200e11 0.73834802771e11 / 0.23471078400e11 0 -0.5274106087e10 / 0.1173553920e10 0.33743985841e11 / 0.5867769600e10 -0.6482602549e10 / 0.2607897600e10 0.1506017269e10 / 0.3911846400e10; -0.76099447e8 / 0.152409600e9 0.63413503537e11 / 0.23471078400e11 -0.1802732209e10 / 0.325987200e9 0.5274106087e10 / 0.1173553920e10 0 -0.7165829063e10 / 0.2607897600e10 0.23903110999e11 / 0.11735539200e11 -0.5346675911e10 / 0.11735539200e11; 0.1397443e7 / 0.12192768e8 -0.4764412871e10 / 0.3911846400e10 0.65514173e8 / 0.16299360e8 -0.33743985841e11 / 0.5867769600e10 0.7165829063e10 / 0.2607897600e10 0 -0.1060918223e10 / 0.11735539200e11 0.628353989e9 / 0.3911846400e10; 0 0.1668252557e10 / 0.5867769600e10 -0.79341409141e11 / 0.58677696000e11 0.6482602549e10 / 0.2607897600e10 -0.23903110999e11 / 0.11735539200e11 0.1060918223e10 / 0.11735539200e11 0 0.25889988599e11 / 0.39118464000e11; 0 -0.842644697e9 / 0.29338848000e11 0.1282384321e10 / 0.7823692800e10 -0.1506017269e10 / 0.3911846400e10 0.5346675911e10 / 0.11735539200e11 -0.628353989e9 / 0.3911846400e10 -0.25889988599e11 / 0.39118464000e11 0;];
 
     Q(1:8,1:8)=Q_U;
-    Q(m-7:m,m-7:m)=flipud( fliplr( -Q_U ) );
+    Q(m-7:m,m-7:m)=rot90(  -Q_U ,2 );
 
-    e_1=zeros(m,1);e_1(1)=1;
-    e_m=zeros(m,1);e_m(m)=1;
+    e_1=sparse(m,1);e_1(1)=1;
+    e_m=sparse(m,1);e_m(m)=1;
 
 
-    D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ;
+    D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -72,13 +76,13 @@
     % M=M/h;
 
     S_U=[-0.12700800e8 / 0.7493827e7 0.185023321e9 / 0.89925924e8 0.39791489e8 / 0.44962962e8 -0.38243253e8 / 0.14987654e8 0.76099447e8 / 0.44962962e8 -0.34936075e8 / 0.89925924e8;]/h;
-    S_1=zeros(1,m);
+    S_1=sparse(1,m);
     S_1(1:6)=S_U;
-    S_m=zeros(1,m);
+    S_m=sparse(1,m);
 
     S_m(m-5:m)=fliplr(-S_U);
 
-    %D2=HI*(-M-e_1*S_1+e_m*S_m);
+    %D2=H\(-M-e_1*S_1+e_m*S_m);
 
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -102,37 +106,40 @@
     %
     %
      S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2;
-     S2_1=zeros(1,m);
+     S2_1=sparse(1,m);
      S2_1(1:5)=S2_U;
-     S2_m=zeros(1,m);
+     S2_m=sparse(1,m);
      S2_m(m-4:m)=fliplr(S2_U);
     %
     %
     %
-    % D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ;
+    % D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ;
 
     % Fourth derivative, 0th order accurate at first 6 boundary points (still
     % yield 4th order convergence if stable: for example u_tt=-u_xxxx
 
     m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8;
-    M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
-
+%     M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
+    stencil = [m4,m3,m2,m1,m0,m1,m2,m3,m4];
+    d = (length(stencil)-1)/2;
+    diags = -d:d;
+    M4 = stripeMatrix(stencil, diags, m);
 
 
     M4_U=[0.600485868980522851e18 / 0.274825314114120000e18 -0.1421010223681841e16 / 0.348984525859200e15 0.38908412970187e14 / 0.1586293299360000e16 0.10224077451922837e17 / 0.2243471951952000e16 -0.7577302712815639e16 / 0.1744922629296000e16 0.138091642084013e15 / 0.59351109840000e14 -0.3775041725375197e16 / 0.4486943903904000e16 0.9907210230881393e16 / 0.61072292025360000e17; -0.1421010223681841e16 / 0.348984525859200e15 0.3985852497808703e16 / 0.407903991264000e15 -0.90048788923861e14 / 0.15579666333000e14 -0.4312795866499e13 / 0.997098645312e12 0.4414634708891947e16 / 0.448694390390400e15 -0.886174803100459e15 / 0.99709864531200e14 0.4333e4 / 0.1000e4 -0.13800578064893047e17 / 0.15704303663664000e17; 0.38908412970187e14 / 0.1586293299360000e16 -0.90048788923861e14 / 0.15579666333000e14 0.2071682582321887e16 / 0.113306664240000e15 -0.769471337294003e15 / 0.41545776888000e14 0.112191585452033e15 / 0.166183107552000e15 0.7204491902193671e16 / 0.623186653320000e15 -0.24847093554379e14 / 0.3115933266600e13 0.943854037768721e15 / 0.545288321655000e15; 0.10224077451922837e17 / 0.2243471951952000e16 -0.4312795866499e13 / 0.997098645312e12 -0.769471337294003e15 / 0.41545776888000e14 0.3086874339649421e16 / 0.81580798252800e14 -0.396009005312111e15 / 0.16618310755200e14 0.348854811893087e15 / 0.249274661328000e15 0.895954627955053e15 / 0.224347195195200e15 -0.184881685054543e15 / 0.166183107552000e15; -0.7577302712815639e16 / 0.1744922629296000e16 0.4414634708891947e16 / 0.448694390390400e15 0.112191585452033e15 / 0.166183107552000e15 -0.396009005312111e15 / 0.16618310755200e14 0.3774861828677557e16 / 0.112173597597600e15 -0.5693689108983593e16 / 0.249274661328000e15 0.803944126167107e15 / 0.99709864531200e14 -0.19547569411550791e17 / 0.15704303663664000e17; 0.138091642084013e15 / 0.59351109840000e14 -0.886174803100459e15 / 0.99709864531200e14 0.7204491902193671e16 / 0.623186653320000e15 0.348854811893087e15 / 0.249274661328000e15 -0.5693689108983593e16 / 0.249274661328000e15 0.73965842628398389e17 / 0.2492746613280000e16 -0.2184472662036043e16 / 0.124637330664000e15 0.46667e5 / 0.10000e5; -0.3775041725375197e16 / 0.4486943903904000e16 0.4333e4 / 0.1000e4 -0.24847093554379e14 / 0.3115933266600e13 0.895954627955053e15 / 0.224347195195200e15 0.803944126167107e15 / 0.99709864531200e14 -0.2184472662036043e16 / 0.124637330664000e15 0.37593640125444199e17 / 0.2243471951952000e16 -0.37e2 / 0.4e1; 0.9907210230881393e16 / 0.61072292025360000e17 -0.13800578064893047e17 / 0.15704303663664000e17 0.943854037768721e15 / 0.545288321655000e15 -0.184881685054543e15 / 0.166183107552000e15 -0.19547569411550791e17 / 0.15704303663664000e17 0.46667e5 / 0.10000e5 -0.37e2 / 0.4e1 0.12766926490502478779e20 / 0.1099301256456480000e19;];
 
     M4(1:8,1:8)=M4_U;
 
-    M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) );
+    M4(m-7:m,m-7:m)=rot90(  M4_U ,2 );
     M4=M4/h^3;
 
     S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3;
-    S3_1=zeros(1,m);
+    S3_1=sparse(1,m);
     S3_1(1:5)=S3_U;
-    S3_m=zeros(1,m);
+    S3_m=sparse(1,m);
     S3_m(m-4:m)=fliplr(-S3_U);
 
-    D4=HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
+    D4=H\(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
 
     S_1 = S_1';
     S_m = S_m';