% below implementation is for S w/ mid-point subdivision
n=20    
tmpDiagonal = zeros(n,1);tmpDiagonal(1)=1;tmpDiagonal(2)=0.5;
    S = spdiags(tmpDiagonal,0,2*n,n); 
    for i = 1:(2*n-1)
        tmpDiagonal = zeros(n,1);tmpDiagonal(i)=0.5;tmpDiagonal(i+1)=1;tmpDiagonal(i+2)=0.5;
        tmpS = spdiags([tmpDiagonal],[-i],2*n,n); 
        S = S + tmpS;    
    end
    tmpDiagonal = zeros(n,1);tmpDiagonal(1)=0.5;
    tmpS = spdiags([tmpDiagonal],[-39],2*n,n);
    S = S + tmpS;