--- a/+sbp/+implementations/d4_compatible_halfvariable_4.m	Fri Sep 09 11:03:13 2016 +0200
+++ b/+sbp/+implementations/d4_compatible_halfvariable_4.m	Fri Sep 09 14:53:41 2016 +0200
@@ -13,45 +13,48 @@
     %h=1/(m-1);
     %h=1;
     
-    BP = 4;
+    BP = 6;
     if(m<2*BP)
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
     end
 
-    c=ones(m,1);
+%     c=ones(m,1);
 
 
-    H=diag(ones(m,1),0);
+    H=speye(m,m);
     H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]);
-    H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48])));
+    H(m-3:m,m-3:m)=rot90(diag([17/48 59/48 43/48 49/48]),2);
     H=H*h;
     HI=inv(H);
     HI = sparse(HI);
 
 
 
-    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.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;];
-    Q(1:4,1:4)=Q_U;
-    Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) );
+%     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));
+      e=ones(m,1);
+%       Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
+%     Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;];
+%     Q(1:4,1:4)=Q_U;
+%     Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,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=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
 
     M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;];
-    M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
+%     M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
+    M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12;
 
     M(1:4,1:4)=M_U;
 
-    M(m-3:m,m-3:m)=flipud( fliplr( M_U ) );
+    M(m-3:m,m-3:m)=rot90(  M_U ,2 );
     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);
     S_1 = S_1';
     S_m = S_m';
@@ -124,15 +127,19 @@
 
 
     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);
     S2_1 = S2_1';
     S2_m = S2_m';
 
     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));
 
@@ -140,13 +147,13 @@
 
     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);
     S3_1 = S3_1';
     S3_m = S3_m';