--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/spdiagsPeriodic.m	Thu Jul 26 18:04:29 2018 -0700
@@ -0,0 +1,60 @@
+function A = spdiagsPeriodic(vals,diags)
+    % Creates an m x m periodic discretization matrix.
+    % vals - m x ndiags matrix of values
+    % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on
+    % vals that are not on main diagonal are going to spill over to
+    % off-diagonal corners.
+
+    default_arg('diags',0);
+
+    [m, ~] = size(vals);
+
+    A = sparse(m,m);
+
+    for i = 1:length(diags)
+
+        d = diags(i);
+        a = vals(:,i);
+
+        % Sub-diagonals
+        if d < 0
+            a_bulk = a(1+abs(d):end);
+            a_corner = a(1:1+abs(d)-1);
+            corner_diag = m-abs(d);
+            A = A + spdiagVariable(a_bulk, d);
+            A = A + spdiagVariable(a_corner, corner_diag);
+
+        % Super-diagonals
+        elseif d > 0
+            a_bulk = a(1:end-d);
+            a_corner = a(end-d+1:end);
+            corner_diag = -m + d;
+            A = A + spdiagVariable(a_bulk, d);
+            A = A + spdiagVariable(a_corner, corner_diag);
+
+        % Main diagonal
+        else
+             A = A + spdiagVariable(a, 0);
+        end
+
+    end
+
+end
+
+function A = spdiagVariable(a,i)
+    default_arg('i',0);
+
+    if isrow(a)
+        a = a';
+    end
+
+    n = length(a)+abs(i);
+
+    if i > 0
+        a = [sparse(i,1); a];
+    elseif i < 0
+        a = [a; sparse(abs(i),1)];
+    end
+
+    A = spdiags(a,i,n,n);
+end