diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h
index 05353f130..e36fbb1d6 100644
--- a/include/openmc/interpolate.h
+++ b/include/openmc/interpolate.h
@@ -4,6 +4,9 @@
 #include <cmath>
 #include <vector>
 
+#include <gsl/gsl-lite.hpp>
+
+#include "openmc/error.h"
 #include "openmc/search.h"
 
 namespace openmc {
@@ -47,7 +50,7 @@ inline double interpolate_lagrangian(gsl::span<const double> xs,
       numerator *= (x - xs[idx + j]);
       denominator *= (xs[idx + i] - xs[idx + j]);
     }
-    output += (numerator / denominator) * ys[i];
+    output += (numerator / denominator) * ys[idx + i];
   }
 
   return output;
diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt
index 50ae9a208..24ec1b7a9 100644
--- a/tests/cpp_unit_tests/CMakeLists.txt
+++ b/tests/cpp_unit_tests/CMakeLists.txt
@@ -2,6 +2,7 @@ set(TEST_NAMES
   test_distribution
   test_file_utils
   test_tally
+  test_interpolate
   # Add additional unit test files here
 )
 
diff --git a/tests/cpp_unit_tests/test_interpolate.cpp b/tests/cpp_unit_tests/test_interpolate.cpp
new file mode 100644
index 000000000..4f19f2b63
--- /dev/null
+++ b/tests/cpp_unit_tests/test_interpolate.cpp
@@ -0,0 +1,51 @@
+#include <iostream>
+#include <map>
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include "openmc/interpolate.h"
+#include "openmc/search.h"
+
+using namespace openmc;
+
+TEST_CASE("Test Lagranian Interpolation")
+{
+  std::vector<double> xs {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+  std::vector<double> ys {0.0, 1.0, 1.0, 2.0, 3.0, 3.0, 5.0};
+
+  // ensure we get data points back at the x values
+  for (int n = 1; n <= 6; n++) {
+    for (int i = 0; i < xs.size(); i++) {
+      double x = xs[i];
+      double y = ys[i];
+
+      size_t idx = lower_bound_index(xs.begin(), xs.end(), x);
+      idx = std::min(idx, xs.size() - n - 1);
+      double out = interpolate_lagrangian(xs, ys, idx, x, n);
+      REQUIRE(out == y);
+    }
+  }
+
+  // spot checks based on an independent implementation of Lagrangian
+  // interpolation
+  std::map<int, std::vector<std::pair<double, double>>> checks;
+  checks[1] = {{0.5, 0.5}, {4.5, 3.0}, {2.5, 1.5}, {5.5, 4.0}};
+  checks[2] = {{2.5, 1.5}, {4.5, 2.75}, {4.9999, 3.0}, {4.00001, 3.0}};
+  checks[3] = {{2.5592, 1.5}, {4.5, 2.9375}, {4.9999, 3.0}, {4.00001, 3.0}};
+
+  for (auto check_set : checks) {
+    int order = check_set.first;
+    auto checks = check_set.second;
+
+    for (auto check : checks) {
+      double input = check.first;
+      double exp_output = check.second;
+
+      size_t idx = lower_bound_index(xs.begin(), xs.end(), input);
+      idx = std::min(idx, xs.size() - order - 1);
+      double out = interpolate_lagrangian(xs, ys, idx, input, order);
+      REQUIRE_THAT(out, Catch::Matchers::WithinAbs(exp_output, 1e-04));
+    }
+  }
+}
\ No newline at end of file