--- a/+sbp/+implementations/d4_compatible_4.m	Fri Sep 09 11:03:13 2016 +0200
+++ b/+sbp/+implementations/d4_compatible_4.m	Fri Sep 09 14:53:41 2016 +0200
@@ -31,10 +31,10 @@
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
     end
 
-    H=diag(ones(m,1),0);
+    H=speye(m,m);
     H_U=[0.3e1 / 0.11e2 0 0 0 0 0; 0 0.2125516311e10 / 0.1311004640e10 0 0 0 0; 0 0 0.278735189e9 / 0.1966506960e10 0 0 0; 0 0 0 0.285925927e9 / 0.163875580e9 0 0; 0 0 0 0 0.1284335339e10 / 0.1966506960e10 0; 0 0 0 0 0 0.4194024163e10 / 0.3933013920e10;];
     H(1:6,1:6)=H_U;
-    H(m-5:m,m-5:m)=fliplr(flipud(H_U));
+    H(m-5:m,m-5:m)=rot90(H_U,2);
     H=H*h;
     HI=inv(H);
 
@@ -42,19 +42,23 @@
     % First derivative SBP operator, 1st order accurate at first 6 boundary points
 
     q2=-1/12;q1=8/12;
-    Q=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=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 = [-q2,-q1,0,q1,q2];
+    d = (length(stencil)-1)/2;
+    diags = -d:d;
+    Q = stripeMatrix(stencil, diags, m);
 
     %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
 
     Q_U = [0 0.9e1 / 0.11e2 -0.9e1 / 0.22e2 0.1e1 / 0.11e2 0 0; -0.9e1 / 0.11e2 0 0.2595224893e10 / 0.2622009280e10 -0.151435707e9 / 0.327751160e9 0.1112665611e10 / 0.2622009280e10 -0.1290899e7 / 0.9639740e7; 0.9e1 / 0.22e2 -0.2595224893e10 / 0.2622009280e10 0 0.1468436423e10 / 0.983253480e9 -0.1194603401e10 / 0.983253480e9 0.72033031e8 / 0.238364480e9; -0.1e1 / 0.11e2 0.151435707e9 / 0.327751160e9 -0.1468436423e10 / 0.983253480e9 0 0.439819541e9 / 0.327751160e9 -0.215942641e9 / 0.983253480e9; 0 -0.1112665611e10 / 0.2622009280e10 0.1194603401e10 / 0.983253480e9 -0.439819541e9 / 0.327751160e9 0 0.1664113643e10 / 0.2622009280e10; 0 0.1290899e7 / 0.9639740e7 -0.72033031e8 / 0.238364480e9 0.215942641e9 / 0.983253480e9 -0.1664113643e10 / 0.2622009280e10 0;];
     Q(1:6,1:6)=Q_U;
-    Q(m-5:m,m-5:m)=flipud( fliplr( -Q_U ) );
+    Q(m-5:m,m-5: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.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
-     S_1=zeros(1,m);
+     S_1=sparse(1,m);
      S_1(1:4)=S_U;
-     S_m=zeros(1,m);
+     S_m=sparse(1,m);
 
      S_m(m-3: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,20 +106,25 @@
 
 
     S2_U=[2 -5 4 -1;]/h^2;
-    S2_1=zeros(1,m);
+    S2_1=sparse(1,m);
     S2_1(1:4)=S2_U;
-    S2_m=zeros(1,m);
+    S2_m=sparse(1,m);
     S2_m(m-3: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
 
     m3=-1/6;m2=2;m1=-13/2;m0=28/3;
-    M4=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=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 = [m3,m2,m1,m0,m1,m2,m3];
+    d = (length(stencil)-1)/2;
+    diags = -d:d;
+    M4 = stripeMatrix(stencil, diags, m);
 
     %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0));
 
@@ -123,16 +132,16 @@
 
     M4(1:6,1:6)=M4_U;
 
-    M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) );
+    M4(m-5:m,m-5:m)=rot90(  M4_U ,2 );
     M4=M4/h^3;
 
     S3_U=[-1 3 -3 1;]/h^3;
-    S3_1=zeros(1,m);
+    S3_1=sparse(1,m);
     S3_1(1:4)=S3_U;
-    S3_m=zeros(1,m);
+    S3_m=sparse(1,m);
     S3_m(m-3: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';