diff --git a/LagrangePolynomial.m b/LagrangePolynomial.m
new file mode 100644
index 0000000..2331a91
--- /dev/null
+++ b/LagrangePolynomial.m
@@ -0,0 +1,29 @@
+function [Lnk] = LagrangePolynomial(x_data, k, x)
+%LagrangePolynomial evaluates a single Lagrange polynomial
+%
+%   Input
+%       - x_data: One-dimensional list of x-coordinates of poins through
+%       which the interpolating polynomial will pass. The degree of the
+%       interpolating polynomial, n, will be inferred from the number
+%       entries in x_data
+%       - k: The interpolating polynomial will equal one at x_data(k + 1)
+%       and be zero for other entries in x_data. 0 <= k <= n
+%        - x: One-dimensional list of values for which the Lagrange
+%        polynomial will be evaluted
+%   Output
+%       - Lnk: Value of the Lagrange polynomial of degree n, index k, at x
+
+n = length(x_data) - 1;
+
+%Make x a row vector. As noted above, it is assumed to be a single row or
+%column
+x = x(:)';
+
+Lnk = 1;
+for(m = 0:n)
+    if(m ~= k)
+        Lnk = Lnk.*(x - x_data(m + 1))/(x_data(k + 1) - x_data(m +1));
+    end
+end
+end
+
diff --git a/TestLagrangePolynomial.m b/TestLagrangePolynomial.m
new file mode 100644
index 0000000..494c575
--- /dev/null
+++ b/TestLagrangePolynomial.m
@@ -0,0 +1,22 @@
+%Test LagrangePolynomial function
+
+close all
+clear all
+
+number_points = 10;
+rng('default')  %Reset the pseudo-random number generator
+
+for(m = 2:number_points)
+    x_data = 10 * rand(1, m);
+    n = m - 1;
+    disp(['Degree = ' num2str(n)])
+    A = zeros(m);
+    for(k = 0:(m - 1))
+        A(k + 1, :) = LagrangePolynomial(x_data, k, x_data);
+    end
+    disp(A)
+    
+    if(~isequal(A, eye(m)))
+        error('Values do not match')
+    end
+end
\ No newline at end of file