diff --git a/src/matlab/alt_pascal.m b/src/matlab/alt_pascal.m
new file mode 100644
index 00000000..8873b975
--- /dev/null
+++ b/src/matlab/alt_pascal.m
@@ -0,0 +1,13 @@
+function [ A ] = alt_pascal( n )
+%   Calculates the alternating pascal triangle up to n rows.
+%   This is the change of basis matrix from the basis {i in {0,1,...,n-1}:(e^(h*D)-1)^i} to 
+%   the standard basis. {i in {0,1,...,n-1}:(e^(i*h*D))}
+    A=eye(n);
+    A(:,1)=2*rem(1:n,2) - 1; %Alternating vector [1,-1,1,-1,...]
+    for i=2:n
+        A(i,2:end)=A(i-1,1:end-1)-A(i-1,2:end);
+    end
+    
+
+
+end
\ No newline at end of file
diff --git a/src/matlab/basisMatrixForGradient.m b/src/matlab/basisMatrixForGradient.m
new file mode 100644
index 00000000..a2cfcd2b
--- /dev/null
+++ b/src/matlab/basisMatrixForGradient.m
@@ -0,0 +1,43 @@
+function [A]=basisMatrixForGradient(k)
+% Returns a k+1 by k+1 triangular matrix.
+% This is the change of basis matrix from the basis {1}U{f_0(x),f_1(x),...,f_k(x)} 
+% to {1,x,x^3,...,x^(2k+1)}.Where f_i(x)=∫_1^x (t^2-1)^idt. 
+% This is the corresponding change of basis matrix for
+% extensionOperatorGradient like altPascal was for extensionOperatorDivergence.
+
+    A=zeros(k+1);
+    A(1,1)=1;
+    currentRow=zeros([1,k+1]);
+    %This is f_i(x), 
+    %Initial row set to f_0(x) which is just x-1
+    currentRow(1)=-1;
+    currentRow(2)=1;
+    alternatingBinom=zeros([1,k+1]); %this list represents x*(x^2-1)^i, which if expanded out
+    % has coordinate vector [0,(-1)^(i),(-1)^(i-1)binom(i,1),...,1]
+    alternatingBinom(2)=-1;
+    alternatingBinom(3)=1;
+    alternatingBinomCopy=zeros([1,k+1]);
+    i=1;
+    leftBoundaryBinom=1; %(-1)^(i), the first term of our alternatingBinom
+    while 1
+        A(i+1,:)=currentRow;
+        if (i==k)
+            break
+        end
+        currentRow=(1/(2*i+1)*(alternatingBinom+currentRow)-currentRow);
+        %f_i(x)=(x(x^2-1)^i)/(2i+1)-(2i)/(2i+1)f_{i-1}(x)
+        alternatingBinomCopy(2)=leftBoundaryBinom;
+        for j=3:(i+2)
+            alternatingBinomCopy(j)=alternatingBinom(j-1)-alternatingBinom(j);
+        end
+        alternatingBinomCopy(i+3)=1;
+        for j=2:(i+3)
+            alternatingBinom(j)=alternatingBinomCopy(j);
+        end
+        leftBoundaryBinom=-leftBoundaryBinom;
+        i=i+1;
+    end
+    
+
+
+end
\ No newline at end of file
diff --git a/src/matlab/calculateInteriorRow.m b/src/matlab/calculateInteriorRow.m
new file mode 100644
index 00000000..eca6384e
--- /dev/null
+++ b/src/matlab/calculateInteriorRow.m
@@ -0,0 +1,45 @@
+function output=calculateInteriorRow(k)
+% Returns the interior row for the kth order mimetic divergence or gradient.
+%
+% Parameters:
+%                k : Order of accuracy
+    output=zeros([1,k]); 
+    % Represents ∑_{i=0}^{k/2-1}
+    % (∏_{j=1}^{i}(0.125/j-0.25))/(2i+1))(e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+1).
+    % Our answer. 
+    r=idivide(k,int32(2));
+    output(r)=-1.0;
+    output(r+1)=1.0; %our starting value for this series is (e^(h/2*d/dx)-e^(-h/2*d/dx))
+    if (k>2)
+        nextTerm=zeros([1,k]);% NextTerm in the series.
+        %Initially set to (e^(h/2*d/dx)-e^(-h/2*d/dx))^(3)
+        nextTermCopy=zeros([1,k]); 
+        nextTerm(r-1)=-1.0;
+        nextTerm(r)=3.0;
+        nextTerm(r+1)=-3.0;
+        nextTerm(r+2)=1.0;
+        scalefactor=1; %this number is ∏_{j=1}^{i}(0.125/j-0.25).
+        i=1;
+        while 1
+            scalefactor=scalefactor*(0.125/i-0.25);
+            output=output+scalefactor/(2.0*i+1.0)*nextTerm;
+            if (i==r-1)
+                break
+            end
+            nextTermCopy(r-i-1)=-1.0; %Updates next term from (e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+1)
+             %to (e^(h/2*d/dx)-e^(-h/2*d/dx))^(2i+3)
+            nextTermCopy(r-i)=1.0*(2.0*i+3.0);
+            nextTermCopy(r+i+1)=-1.0*(2.0*i+3.0);
+            nextTermCopy(r+i+2)=1.0;
+            for j=(r-i+1):(r+i)
+                nextTermCopy(j)=nextTerm(j+1)-2*nextTerm(j)+nextTerm(j-1);
+            end
+            for j=(r-i-1):(r+i+2)
+                nextTerm(j)=nextTermCopy(j);
+            end
+            i=i+1;
+
+        
+        end
+    end
+end
\ No newline at end of file
diff --git a/src/matlab/div.m b/src/matlab/div.m
index 9c1e95b5..6f8e8779 100644
--- a/src/matlab/div.m
+++ b/src/matlab/div.m
@@ -26,12 +26,14 @@
     assert(k >= 2, 'k >= 2');
     assert(mod(k, 2) == 0, 'k % 2 = 0');
     assert(m >= 2*k+1, ['m >= ' num2str(2*k+1) ' for k = ' num2str(k)]);
+    n_rows=m+2;
+    n_cols=m+1;
 
-    D = sparse(m+2, m+1);
+    D = sparse(n_rows,n_cols);
 
     switch k
         case 2
-            for i = 2:m+1
+            for i = 2:n_cols
                D(i, i-1:i) = [-1 1];
             end
 
@@ -39,7 +41,7 @@
             A = [-11/12 17/24 3/8 -5/24 1/24];
             D(2, 1:5) = A;
             D(m+1, m-3:end) = -fliplr(A);
-            for i = 3:m
+            for i = 3:n_cols-1
                D(i, i-2:i+1) = [1/24 -9/8 9/8 -1/24];
             end
 
@@ -48,7 +50,7 @@
                     31/960  -687/640 129/128   19/192 -3/32    21/640  -3/640];
             D(2:3, 1:7) = A;
             D(m:m+1, m-5:end) = -rot90(A,2);
-            for i = 4:m-1
+            for i = 4:n_cols-2
                 D(i, i-3:i+2) = [-3/640 25/384 -75/64 75/64 -25/384 3/640];
             end
 
@@ -58,9 +60,60 @@
                    -59/17920    1175/21504 -1165/1024   1135/1024    25/3072  -251/5120   25/1024   -45/7168     5/7168];
             D(2:4, 1:9) = A;
             D(m-1:m+1, m-7:end) = -rot90(A,2);
-            for i = 5:m-2
+            for i = 5:n_cols-3
                 D(i, i-4:i+3) = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168];
             end
+        otherwise 
+            neighbors = zeros(1, k); % Bandwidth = k
+            neighbors(1) = 1/2 - k/2;
+            for i = 2 : k
+                neighbors(i) = neighbors(i-1)+1;
+            end
+    
+    
+            A = invvander(neighbors);
+    
+   
+            b = zeros(k, 1);
+            b(2) = 1;
+    
+    
+            coeffs = A*b;
+    
+            j = 1;
+            for i = k/2+1 : n_rows - k/2
+                D(i, j:j+k-1) = coeffs;
+                j = j + 1;
+            end
+    
+            p = k/2-1;
+            q = k+1;
+            A = sparse(p, q);
+            for i = 1 : p % For each row of A
+                neighbors = zeros(1, q); % k+1 points are used for the boundaries
+                neighbors(1) = 1/2 - i; % Shifting the stencil to the right
+                for j = 2 : q
+                    neighbors(j) = neighbors(j-1)+1;
+                end
+                V = invvander(neighbors);
+                b = zeros(q, 1);
+                b(2) = 1;
+                coeffs = V*b;
+                A(i, 1:q) = coeffs;
+            end
+   
+            D(2:p+1, 1:q) = A;
+    
+    
+
+            Pp = fliplr(speye(p));
+            Pq = fliplr(speye(q));
+        
+            A = -Pp*A*Pq;
+    
+    
+           D(n_rows-p:n_rows-1, n_cols-q+1:n_cols) = A;
     end
     D = (1/dx).*D;
+
 end
diff --git a/src/matlab/div1dfft.m b/src/matlab/div1dfft.m
new file mode 100644
index 00000000..5412bda3
--- /dev/null
+++ b/src/matlab/div1dfft.m
@@ -0,0 +1,11 @@
+function [du]=div1dfft(u,k,dx)
+    m=size(u)
+    m=m(1)
+    A=extensionOperatorDivergence(k,m-1);
+    deriv=calculateInteriorRow(k)'; %size is k-1
+  
+    uExt=A*u; %size is m+k-2
+    du=-ifft(fft(deriv,m+2*k-4) .* fft(uExt,m+2*k-4))/dx;
+    du=du(k:(k+m-2))
+    du=[0;du;0]
+end
\ No newline at end of file
diff --git a/src/matlab/divergencehelper.m b/src/matlab/divergencehelper.m
new file mode 100644
index 00000000..a4c23ecd
--- /dev/null
+++ b/src/matlab/divergencehelper.m
@@ -0,0 +1,25 @@
+function D=divergencehelper(k,m)
+%Calculates the mimetic divergence operator of order k,
+%   barring dividing by dx. 
+%   k :order of accuracy.
+%   m :number of cells. 
+    interior=calculateInteriorRow(k); %The interior row of our matrix 
+    numNonZeros=m*k; % The number of nonZeroElements in the interior stencil matrix
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    for i=2:(m+1)
+        for j=1:(k)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=i+j-2;
+            valueList(elementsInserted)=interior(j);
+            elementsInserted=elementsInserted+1;
+        end
+    end
+    interiorScheme=sparse(rowList,columnList,valueList,m+2,m+k-1);
+    % This matrix is just using the interior scheme.
+    % Since we've extended our vector field beyond the boundaries it'd normally have, we don't need to use
+    % seperate power series expansions for the stencils near the boundary. 
+    D=interiorScheme*extensionOperatorDivergence(k,m);
+end
\ No newline at end of file
diff --git a/src/matlab/extensionOperatorDivergence.m b/src/matlab/extensionOperatorDivergence.m
new file mode 100644
index 00000000..77cd69e1
--- /dev/null
+++ b/src/matlab/extensionOperatorDivergence.m
@@ -0,0 +1,53 @@
+function E=extensionOperatorDivergence(k,m)
+% Returns a m+k-1 by m+1 one matrix. Which pads an m+1` point vector field
+% With k/2-1 additional points on each side. These points are
+% approximations of
+% f(-(k/2-1)h),f(-(k/2-2)h),....,f(-h),f(1+h),f(1+2h),...,f(1+(k/2-1)h),
+% with O(h^(k+1)) error on each of the approximations. 
+% 
+% 
+% Parameters:
+%                k : Order of accuracy
+%                m : Number of cells
+    numOfRows=m+k-1;
+    numOfColumns=m+1;
+    numNonZeros=(k-2)*(k+1)+m+1;%The number of nonzero elements in E.
+    rowList=zeros([1,numNonZeros]);%Row indices for non zero elements
+    columnList=zeros([1,numNonZeros]); %Column indices for non zero elements
+    valueList=zeros([1,numNonZeros]); %Values of non zero elements
+    elementsInserted=1; %Number of elements inserted into rowList,ColumnList
+    %and ValueList
+    r=idivide(k,int32(2));
+    R=zeros(r-1,k+1);% Calculates the matrix R_{i,j}=binom(i-k/2,j-1) 
+    R(:,1)=ones([r-1,1]);
+    currentTerm=(-1);
+    for i=2:(k+1)
+        R(r-1,i)=currentTerm;
+        currentTerm=currentTerm*(-1);
+    end
+    for i=(r-2):(-1):1
+        for j=2:(k+1)
+            R(i,j)=R(i+1,j)-R(i,j-1);
+        end
+    end
+    R=R*alt_pascal(k+1); %Converts R into the standard basis 
+    for i=1:(r-1) %
+        for j=1:(k+1)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=j;
+            valueList(elementsInserted)=R(i,j);
+            rowList(elementsInserted+1)=m+k-i;
+            columnList(elementsInserted+1)=m+2-j;
+            valueList(elementsInserted+1)=R(i,j);
+            elementsInserted=elementsInserted+2;
+        end
+    end
+    for i=r:m+r
+        rowList(elementsInserted)=i;
+        columnList(elementsInserted)=i+1-r;
+        valueList(elementsInserted)=1;
+        elementsInserted=elementsInserted+1;
+    end
+    E=sparse(rowList,columnList,valueList,numOfRows,numOfColumns);
+
+end
\ No newline at end of file
diff --git a/src/matlab/extensionOperatorGradient.m b/src/matlab/extensionOperatorGradient.m
new file mode 100644
index 00000000..e57401fc
--- /dev/null
+++ b/src/matlab/extensionOperatorGradient.m
@@ -0,0 +1,48 @@
+function E=extensionOperatorGradient(k,m)
+    % Returns a m+k by m+2 matrix, 
+    % Which turns the list [f(0),f(h/2),f(3h/2),...,f(1-h/2),f(1)]
+    % Into
+    % [f(-(k/2-1/2)h),f(-(k/2-3/2)h),...,f(-h/2),
+    % f(h/2),...f(1-h/2),f(1+h/2),...f(1+(k/2-1/2)*h)]
+    numOfRows=m+k;
+    numOfColumns=m+2;
+    r=idivide(k,int32(2));
+    numNonZeros=m+k*(k+1); 
+    % The number of nonZeroElements in the interior stencil matrix
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    R=zeros([r,k+1]);
+    R(:,2)=ones([1,r]);
+    R(r,2:end)=2*rem(1:k,2) - 1;
+    for i=(r-1):(-1):1
+        for j=3:(k+1)
+            R(i,j)=R(i+1,j)-R(i,j-1);
+        end
+    end
+
+    Scalefactor=diag(-(k-1):2:-1);
+    R=Scalefactor*R;
+    R(:,1)=ones([1,r]);
+    R=R*basisMatrixForGradient(k);
+    for i=1:r
+        for j=1:(k+1)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=j;
+            valueList(elementsInserted)=R(i,j);
+            rowList(elementsInserted+1)=m+k+1-i;
+            columnList(elementsInserted+1)=m+3-j;
+            valueList(elementsInserted+1)=R(i,j);
+            elementsInserted=elementsInserted+2;
+        end
+    end
+    for j=2:(m+1)
+        rowList(elementsInserted)=r+j-1;
+        columnList(elementsInserted)=j;
+        valueList(elementsInserted)=1;
+        elementsInserted=elementsInserted+1;
+    end
+    E=sparse(rowList,columnList,valueList,numOfRows,numOfColumns);
+
+end
\ No newline at end of file
diff --git a/src/matlab/grad.m b/src/matlab/grad.m
index 6d5edf54..aa270d29 100644
--- a/src/matlab/grad.m
+++ b/src/matlab/grad.m
@@ -67,6 +67,8 @@
             for i = 5:m-3
                 G(i, i-3:i+4) = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168];
             end
+        otherwise 
+            G=gradienthelper(k,m);
     end
     G = (1/dx).*G;
 end
diff --git a/src/matlab/gradienthelper.m b/src/matlab/gradienthelper.m
new file mode 100644
index 00000000..6e5e5af2
--- /dev/null
+++ b/src/matlab/gradienthelper.m
@@ -0,0 +1,19 @@
+function G=gradienthelper(k,m)
+    interior=calculateInteriorRow(k);
+    numNonZeros=(m+1)*k;
+    rowList=zeros([1,numNonZeros]);
+    columnList=zeros([1,numNonZeros]);
+    valueList=zeros([1,numNonZeros]);
+    elementsInserted=1;
+    for i=1:(m+1)
+        for j=1:(k)
+            rowList(elementsInserted)=i;
+            columnList(elementsInserted)=i+j-1;
+            valueList(elementsInserted)=interior(j);
+            elementsInserted=elementsInserted+1;
+        end
+    end
+    interiorScheme=sparse(rowList,columnList,valueList);
+    G=interiorScheme*extensionOperatorGradient(k,m);
+
+end
\ No newline at end of file
diff --git a/src/matlab/invvander.m b/src/matlab/invvander.m
new file mode 100644
index 00000000..d8d34bcf
--- /dev/null
+++ b/src/matlab/invvander.m
@@ -0,0 +1,56 @@
+function B = invvander(v, m)
+%INVVANDER inverse or pseudoinverse of transpose of vandermode matrix 
+
+
+assert(isrow(v), 'v must be a row vector.')
+assert(numel(v) == numel(unique(v)), 'all elements in the v must be distinct.')
+n = numel(v);
+if nargin == 1
+    m = n;
+else
+    assert(isscalar(m), 'm must be a scalar.')
+    assert(m > 0 && mod(m, 1) == 0, 'm must be a poistive integer.')
+end
+
+if m == n
+    % 1-by-1 matrix
+    if n == 1
+        B = 1 / v;
+        return;
+    end
+    
+    p = [-v(1) 1];
+    C(1) = 1;
+    for i = 2:n
+        p = [0, p] + [-v(i) * p, 0];
+        
+        Cp = C;
+        C = zeros(1, i);
+        for j = 1:i - 1
+            C(j) = Cp(j) / (v(j) - v(i));       
+        end
+        C(i) = - sum(C);
+    end
+    
+    B = zeros(n);
+    c = zeros(1, n);
+    
+    for i = 1:n
+        c(n) = 1;
+        for j = n-1:-1:1
+            c(j) = p(j + 1) + v(i) * c(j + 1);
+        end
+        B(i, :) = c * C(i);
+    end
+    
+else
+    V = (v.' .^ (0:(m - 1)));
+    if m > n % over-determined
+        V = V.';
+    end
+    [~, R] = qr(V);
+    B = mldivide(R, mldivide(R', V'));
+    if m < n % under-determined
+        B = B.';
+    end
+end
\ No newline at end of file
diff --git a/tests/matlab/run_tests.m b/tests/matlab/run_tests.m
index 823fee3d..fa5ad170 100644
--- a/tests/matlab/run_tests.m
+++ b/tests/matlab/run_tests.m
@@ -21,3 +21,9 @@
 
 disp('Running: BC Ops consistency test...');
 run('test6.m');
+
+disp('Running: Correctness test of gradienthelper...');
+run('test7.m');
+
+disp('Running: Correctness test of divergencehelper');
+run('test8.m');
\ No newline at end of file
diff --git a/tests/matlab/test7.m b/tests/matlab/test7.m
new file mode 100644
index 00000000..065df19c
--- /dev/null
+++ b/tests/matlab/test7.m
@@ -0,0 +1,23 @@
+% Correctness test of divergence
+
+addpath('../../src/matlab')
+
+tol = 1e-4;
+
+testfun=@ (x) (x-x^2)^2;
+deriv=@ (x) (2*(1-2*x)*(x-x^2));
+points=10;
+for i=10:2:(10+(2*points))
+    faces=linspace(0,1,2*i+2)';
+    B=div(i,2*i+1,1/(2*i+1));
+    centers=[0; linspace(0.5/(2*i+1),1-0.5/(2*i+1),2*i+1)' ;1];
+    testoutput=arrayfun(deriv,centers);
+    d=B*(arrayfun(testfun,faces));
+    if (norm(d-testoutput,inf) > tol)
+        fprintf("Test FAILED!\n");
+        fprintf("error between numerical and actual result is\n");
+        norm(d-testoutput,inf)     
+        i
+    end  
+end
+
diff --git a/tests/matlab/test8.m b/tests/matlab/test8.m
new file mode 100644
index 00000000..42b1449d
--- /dev/null
+++ b/tests/matlab/test8.m
@@ -0,0 +1,23 @@
+% Correctness test of gradient
+
+addpath('../../src/matlab')
+
+tol = 1e-4;
+
+testfun=@ (x) (x-x^2)^2;
+deriv=@ (x) (2*(1-2*x)*(x-x^2));
+points=10;
+for i=10:2:(10+(2*points))
+    faces=linspace(0,1,2*i+1)';
+    B=grad(i,2*i,1/(2*i));
+    centers=[0; linspace(0.5/(2*i),1-0.5/(2*i),2*i)' ;1];
+    testoutput=arrayfun(deriv,faces);
+    d=B*(arrayfun(testfun,centers));
+    if (norm(d-testoutput,inf) > tol)
+        fprintf("Test FAILED!\n");
+        fprintf("error between numerical and actual result is\n");
+        norm(d-testoutput,inf)     
+        i;
+    end  
+end
+fprintf("Test PASSED!\n");
\ No newline at end of file