diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py
index ead6a60..4ead61b 100644
--- a/src/mvesuvio/analysis_reduction.py
+++ b/src/mvesuvio/analysis_reduction.py
@@ -12,7 +12,7 @@
                             VesuvioThickness, Integration, Divide, Multiply, DeleteWorkspaces, \
                             CreateWorkspace, CreateSampleWorkspace
 
-from mvesuvio.util.analysis_helpers import loadConstants, numericalThirdDerivative
+from mvesuvio.util.analysis_helpers import loadConstants, numerical_third_derivativ, extend_range_of_array
 
 
 
@@ -661,18 +661,26 @@ def _neutron_compton_profiles(self, pars):
         gaussRes, lorzRes = self.caculateResolutionForEachMass(centers)
         totalGaussWidth = np.sqrt(widths**2 + gaussRes**2)
 
-        JOfY = scipy.special.voigt_profile(self._y_space_arrays[self._row_being_fit] - centers, totalGaussWidth, lorzRes)
+        # Extend range because third derivative trims it
+        y_space_arrays_extended = extend_range_of_array(self._y_space_arrays[self._row_being_fit], 6)
+
+        JOfY = scipy.special.voigt_profile(y_space_arrays_extended - centers, totalGaussWidth, lorzRes)
 
         FSE = (
-            -numericalThirdDerivative(self._y_space_arrays[self._row_being_fit], JOfY)
+            -numerical_third_derivative(y_space_arrays_extended, JOfY)
             * widths**4
             / deltaQ
             * 0.72
         )
-        scaling_factor = intensities * E0 * E0 ** (-0.92) * masses / deltaQ
-        JOfY *= scaling_factor
-        FSE *= scaling_factor
-        return JOfY+FSE, FSE
+        # Trim profile to match FSE shape
+        JOfY = JOfY[:, 6:-6]
+        return intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ
+
+
+    def _extend_range_of_array(self, arr, n_columns):
+        left_extend = arr[:, :n_columns] + (arr[:, 0] - arr[:, n_columns]).reshape(-1, 1)
+        right_extend = arr[:, -n_columns:] + (arr[:, -1] - arr[:, -n_columns-1]).reshape(-1, 1)
+        return np.concatenate([left_extend, arr, right_extend], axis=-1)
 
 
     def caculateResolutionForEachMass(self, centers):
diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py
index fa2c974..cc074d4 100644
--- a/src/mvesuvio/util/analysis_helpers.py
+++ b/src/mvesuvio/util/analysis_helpers.py
@@ -249,7 +249,13 @@ def loadConstants():
     return constants
 
 
-def numericalThirdDerivative(x, y):
+def extend_range_of_array(arr, n_columns):
+    left_extend = arr[:, :n_columns] + (arr[:, 0] - arr[:, n_columns]).reshape(-1, 1)
+    right_extend = arr[:, -n_columns:] + (arr[:, -1] - arr[:, -n_columns-1]).reshape(-1, 1)
+    return np.concatenate([left_extend, arr, right_extend], axis=-1)
+
+
+def numerical_third_derivative(x, y):
     k6 = (- y[:, 12:] + y[:, :-12]) * 1
     k5 = (+ y[:, 11:-1] - y[:, 1:-11]) * 24
     k4 = (- y[:, 10:-2] + y[:, 2:-10]) * 192
@@ -257,13 +263,9 @@ def numericalThirdDerivative(x, y):
     k2 = (+ y[:, 8:-4] - y[:, 4:-8]) * 387
     k1 = (- y[:, 7:-5] + y[:, 5:-7]) * 1584
 
-    dev = k1 + k2 + k3 + k4 + k5 + k6
-    dev /= np.power(x[:, 7:-5] - x[:, 6:-6], 3)
-    dev /= 12**3
-
-    derivative = np.zeros_like(y)
-    # Padded with zeros left and right to return array with same shape
-    derivative[:, 6:-6] = dev
+    derivative = k1 + k2 + k3 + k4 + k5 + k6
+    derivative /= np.power(x[:, 7:-5] - x[:, 6:-6], 3)
+    derivative /= 12**3
     return derivative
 
 
diff --git a/tests/unit/analysis/test_analysis_helpers.py b/tests/unit/analysis/test_analysis_helpers.py
index d6b20e5..3f60e61 100644
--- a/tests/unit/analysis/test_analysis_helpers.py
+++ b/tests/unit/analysis/test_analysis_helpers.py
@@ -1,11 +1,13 @@
 import unittest
 import numpy as np
+import scipy
 import dill
 import numpy.testing as nptest
 from mock import MagicMock
 from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table,  \
-    fix_profile_parameters, calculate_h_ratio
+    fix_profile_parameters, calculate_h_ratio, numerical_third_derivative, extend_range_of_array
 from mantid.simpleapi import CreateWorkspace, DeleteWorkspace
+import matplotlib.pyplot as plt
 
 
 class TestAnalysisHelpers(unittest.TestCase):
@@ -126,5 +128,28 @@ def test_conversion_of_constraints(self):
         self.assertEqual(converted_constraints[1]['fun']([0, 0, 0, 2, 0, 0, 1]), 2-0.7234)
 
 
+    def test_extend_range_of_array_for_increasing_range(self):
+        x = np.arange(10)
+        x = np.vstack([x, 2*x])
+        x_extended = extend_range_of_array(x, 5)
+        np.testing.assert_array_equal(x_extended, np.vstack([np.arange(-5, 15, 1), np.arange(-10, 30, 2)]))
+
+
+    def test_extend_range_of_array_for_decreasing_range(self):
+        x = np.linspace(-5, 5, 21)
+        x = np.vstack([x, 2*x])
+        x_extended = extend_range_of_array(x, 5)
+        np.testing.assert_array_equal(x_extended, np.vstack([np.linspace(-7.5, 7.5, 31), np.linspace(-15, 15, 31)]))
+
+
+    def test_numerical_third_derivative(self):
+        x= np.linspace(-20, 20, 300)    # Workspaces are about 300 points of range
+        x = np.vstack([x, 2*x])
+        y = scipy.special.voigt_profile(x, 5, 5)
+        numerical_derivative = numerical_third_derivative(x, y)
+        expected_derivative = np.array([np.gradient(np.gradient(np.gradient(y_i, x_i), x_i), x_i)[6: -6] for y_i, x_i in zip(y, x) ])
+        np.testing.assert_allclose(numerical_derivative, expected_derivative, atol=1e-6)
+
+
 if __name__ == "__main__":
     unittest.main()