--- a/+sbp/higher_variable4.m	Mon Sep 05 16:34:22 2016 +0200
+++ b/+sbp/higher_variable4.m	Mon Sep 05 16:38:41 2016 +0200
@@ -7,56 +7,61 @@
     %%% D2          (approx andra derivatan)    %%%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+    Hv = ones(m,1);
+    Hv(1:4) = [17/48 59/48 43/48 49/48];
+    Hv(m-3:m) = rot90(Hv(1:4),2);
+    Hv = h * Hv;
+    H = spdiag(Hv,0);
+    HI = spdiag(1./Hv,0);
 
-    %m=20; %problemstorlek
-    %h=1/(m-1);
-    %h=1;
+    e_1 = sparse(m, 1);
+    e_1(1) = 1;
+    e_m = rot90(e_1, 2);
+
+    S_1 = sparse(m, 1);
+    S_1(1:4) = [-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1]/h;
+    S_m = -rot90(S_1, 2);
 
-    c = ones(m,1);
+    S2_1 = sparse(m, 1);
+    S2_1(1:4) = [2 -5 4 -1;]/h^2;
+    S2_m = rot90(S2_1, 2);
+
+    S3_1 = sparse(m, 1);
+    S3_1(1:4) = [-1 3 -3 1;]/h^3;
+    S3_m = rot90(S3_1, 2);
 
 
-    H=diag(ones(m,1),0);
-    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=H*h;
-    HI=inv(H);
-    HI = sparse(HI);
-
+    q_diags = [-2, -1, 0, 1, 2];
+    q_stencil   = [-1/12, 8/12, 0, -8/12, 1/12];
+    Q = stripeMatrix(q_stencil, q_diags,m);
 
+    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/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(m-3:m,m-3:m)=-rot90(Q_U, 2);
 
-    e_1=zeros(m,1);e_1(1)=1;
-    e_m=zeros(m,1);e_m(m)=1;
 
     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_diags = [-2, -1, 0, 1, 2];
+    m_stencil   = [-1/12, 16/12, -30/12, -16/12, -1/12];
+
+    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: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(1:4)=S_U;
-    S_m=zeros(1,m);
-    S_m(m-3:m)=fliplr(-S_U);
-    S_1 = S_1';
-    S_m = S_m';
-
-
-    M=sparse(m,m);
-    e_1 = sparse(e_1);
-    e_m = sparse(e_m);
-    S_1 = sparse(S_1);
-    S_m = sparse(S_m);
-
 
     scheme_width = 5;
     scheme_radius = (scheme_width-1)/2;
@@ -116,15 +121,6 @@
     end
     D2 = @D2_fun;
 
-
-    S2_U=[2 -5 4 -1;]/h^2;
-    S2_1=zeros(1,m);
-    S2_1(1:4)=S2_U;
-    S2_m=zeros(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);
 
@@ -137,13 +133,6 @@
     M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) );
     M4=M4/h^3;
 
-    S3_U=[-1 3 -3 1;]/h^3;
-    S3_1=zeros(1,m);
-    S3_1(1:4)=S3_U;
-    S3_m=zeros(1,m);
-    S3_m(m-3:m)=fliplr(-S3_U);
-    S3_1 = S3_1';
-    S3_m = S3_m';
 
     D4=HI*(M4-e_1*S3_1'+e_m*S3_m'  + S_1*S2_1'-S_m*S2_m');