diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h
index 0683140339..2764b7ffc2 100644
--- a/include/openmc/interpolate.h
+++ b/include/openmc/interpolate.h
@@ -34,41 +34,11 @@ inline double interpolate_log_log(
   return y0 * exp(f * log(y1 / y0));
 }
 
-inline double interpolate(const std::vector<double>& xs,
-  const std::vector<double>& ys, int idx, double x,
-  Interpolation i = Interpolation::lin_lin)
-{
-  switch (i) {
-  case Interpolation::lin_lin:
-    return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
-  case Interpolation::log_log:
-    return interpolate_log_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
-  case Interpolation::lin_log:
-    return interpolate_lin_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
-  case Interpolation::log_lin:
-    return interpolate_log_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
-  default:
-    fatal_error("Unsupported interpolation");
-  }
-}
-
-inline double interpolate(const std::vector<double>& xs,
-  const std::vector<double>& ys, double x,
-  Interpolation i = Interpolation::lin_lin)
-{
-  int idx = lower_bound_index(xs.begin(), xs.end(), x);
-  return interpolate(xs, ys, idx, x, i);
-}
-
 inline double interpolate_lagrangian(
-  const std::vector<double>& xs, const std::vector<double>& ys, double x)
+  const std::vector<double>& xs, const std::vector<double>& ys, int idx, double x, int order)
 {
-  int idx = lower_bound_index(xs.begin(), xs.end(), x);
-
   std::vector<double> coeffs;
 
-  int order = 3;
-
   for (int i = 0; i < order + 1; i++) {
     double numerator {1.0};
     double denominator {1.0};
@@ -85,6 +55,37 @@ inline double interpolate_lagrangian(
     coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0);
 }
 
+inline double interpolate(const std::vector<double>& xs,
+  const std::vector<double>& ys, double x,
+  Interpolation i = Interpolation::lin_lin)
+{
+  int idx = lower_bound_index(xs.begin(), xs.end(), x);
+
+  switch (i) {
+  case Interpolation::lin_lin:
+    return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
+  case Interpolation::log_log:
+    return interpolate_log_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
+  case Interpolation::lin_log:
+    return interpolate_lin_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
+  case Interpolation::log_lin:
+    return interpolate_log_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
+  case Interpolation::quadratic:
+    // move back one point if x is in the last interval of the x-grid
+    if (idx == xs.size() - 2 && idx > 0) idx--;
+    return interpolate_lagrangian(xs, ys, x, idx, 2);
+  case Interpolation::cubic:
+    // if x is not in the first interval of the x-grid, move back one
+    if (idx > 0) idx--;
+    // if the index was the last interval of the x-grid, move it back one more
+    if (idx == xs.size() - 3) idx--;
+    return interpolate_lagrangian(xs, ys, x, idx, 3);
+  default:
+    fatal_error("Unsupported interpolation");
+  }
+}
+
+
 } // namespace openmc
 
 #endif
\ No newline at end of file
diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp
index f883975a36..a09854e847 100644
--- a/src/tallies/filter_energyfunc.cpp
+++ b/src/tallies/filter_energyfunc.cpp
@@ -75,31 +75,8 @@ void EnergyFunctionFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
   if (p.E_last() >= energy_.front() && p.E_last() <= energy_.back()) {
-    // Search for the incoming energy bin.
-    auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last());
 
-    double f, w;
-    switch (interpolation_) {
-    case Interpolation::lin_lin:
-      w = interpolate_lin_lin(
-        energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
-      break;
-    case Interpolation::lin_log:
-      w = interpolate_lin_log(
-        energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
-      break;
-    case Interpolation::log_lin:
-      w = interpolate_log_lin(
-        energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
-      break;
-    case Interpolation::log_log:
-      w = interpolate_log_log(
-        energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
-      break;
-    default:
-      fatal_error(fmt::format(
-        "Invalid interpolation scheme found on EnergyFunctionFilter {}", id()));
-    }
+    double w = interpolate(energy_, y_, p.E_last(), interpolation_);
 
     // Interpolate on the lin-lin grid.
     match.bins_.push_back(0);